|Home | About | Journals | Submit | Contact Us | Français|
Molecular dynamics (MD) has become a popular method to study ion channels by theoretical means and to provide new insights into their fundamental properties, such as fast conduction and ion specificity. This Perspective deals with one of the current challenges of biomolecular MD studies: the accurate description of polarization. Polarization can be defined as the spatial changes in charge distribution due to the presence of an electric field. In the case of ion channels, theoretical studies ought to be able to describe the response of electronic clouds to the moving ions (i.e., ion-induced polarization). However, the most widely used empirical potential energy functions (force fields) such as AMBER, CHARMM, and GROMOS, are not explicitly polarizable, but rather include polarizability implicitly in an average way in their parameterization. Polarizable extensions of these force fields are under active development, but the lack of in situ reference data that can be used to assess the performance of these force fields renders their development more difficult. Here, we discuss the possible use of quantum mechanics (QM)/molecular mechanics (MM) simulations to assist the development of improved force fields for ion channel studies. These QM/MM simulations highlight some of the possible deficiencies of current force fields and provide examples of the importance of polarization for the accurate description of ion conduction and selectivity.
Ion channels are involved in a broad range of the most basic biological processes, from excitation and signaling to secretion and absorption (Hille, 2001), and are important drug targets (e.g., see “Perspectives on How to Drug an Ion Channel”; May 2008; JGP 131:395). Ion channel modulators are used to treat a variety of diseases, such as epilepsy, hypertension, diabetes, and chronic pain. In addition, it was found recently that K channels are overexpressed in some cancer cells, which has prompted interest in K channels as targets of new anticancer drugs (Kunzelmann, 2005). Presently, most of the drugs that are on the market had to be developed without extensive knowledge of the molecular interactions between the drugs and their targets. However, several atomic structures of ion channels have been elucidated in recent years, and there is some hope that in the future, computer-assisted drug design might have a positive impact on this process. A major prerequisite for such an endeavor, however, is an accurate computational procedure to evaluate drug–channel interactions.
In the past decade, the availability of high-resolution structures for a variety of ion channels has lead to a leap in the number of systems that have been studied with computer simulations. MD has become the most popular theoretical method for ion channel studies because it can provide an atomistic picture of ion conduction and help rationalize the mechanisms involved in ion specificity. In addition, the different “kinetic states” of ion channels and the gating process can be studied in high spatial and temporal resolution. A recent review by Khalili-Araghi et al. (2009) discusses simulation studies of ion channels and transporters that have been performed in the past two years.
Classical MD is based on the assumption that interatomic interactions can be described accurately with simple potential energy functions. Despite the many successes of MD simulations in providing an atomic-level picture of biological processes (Jensen et al., 2010), the method does not produce accurate results in all cases. One of the most serious challenges in this respect is the treatment of electronic polarization because most current force fields are based on a fixed-point charge description of the electrostatic field and thus only capture polarization effects in a mean-field way. A prime example of the importance of polarization is the change of the dipole moment of water in different environments. In gas phase, the water dipole is 1.85 D, but in condensed phase, it reaches a value between 2.4 and 2.6 D, as suggested from classical MD simulations of the dielectric properties, whereas ab initio MD (AIMD) simulations and analysis of experimental data (for a discussion see Lopes et al., 2009) suggest an even larger increase in dipole moment to 2.95 D. Nonpolarizable force fields, on the other hand, typically assign a fixed dipole to water molecules, which is slightly higher than the gas phase dipole, but lower than the maximum value for a dipole in water. For example, the popular TIP3 model uses a fixed dipole moment of 2.35 D. Although this approximation is valid in most cases, it may fail for systems where water molecules travel through very different environments. In addition, in ion channels, the fundamental interactions between ions and the protein are essentially electrostatic in nature and cannot be described accurately without electronic polarization.
Concerns about the description of polarization in MD were already raised in the 1970s, and a growing interest for methods that are able to provide a more accurate description of polarization can be noted in the current literature. A recent review has focused specifically on the polarization effects in K channels (Illingworth and Domene, 2009). Several other excellent reviews have been published that discuss polarization effects in MD. Ponder and Case (2003), for example, discuss the development of force fields and provide some historical background. Yu and van Gunsteren (2005) review various ways to implement polarization into MD programs. Warshel et al. (2007) supply useful benchmarks that can be exploited to develop polarizable force fields and discuss examples of biological systems that are not well described within a nonpolarizable approach. Finally, Cieplak et al. (2009) provide an up-to-date overview of the latest developments in this field. Although polarizable extensions of widely used force fields (AMBER, CHARMM, and GROMOS) are underway, their applications to macromolecules of biological importance are still limited, and the methods and software required to treat polarization are not yet as standardized as for nonpolarizable protein potentials.
In recent years, a radically different approach has emerged to study ion channels: the use of ab initio simulations. In AIMD, the interatomic forces are obtained by solving the electronic structure problem in the framework of QM. Density functional theory (DFT) is a successful QM theory that is often used, in conjunction with AIMD, to solve the electronic structure problem on the fly. DFT was not considered accurate enough for calculations in quantum chemistry until the 1990s, when the approximations used in the theory were greatly refined to better model the exchange and correlation interactions. DFT is now a leading method for electronic structure calculations in chemistry, solid-state physics, and the most commonly applied method for quantum mechanical calculations in biological systems because it often provides an excellent balance between computational cost and accuracy. However, despite improvements in the method, there are still difficulties associated with the proper description of some intermolecular interactions with DFT, especially in the case of van der Waals forces, charge transfer excitations, and strongly correlated systems. Functionals can be improved by including a dependence on density gradients, as in the generalized gradient approximation (GGA), and by including kinetic energy density terms (meta-GGA approaches). Although meta-GGAs and hybrid meta–GGAs are frequently more accurate on average than GGAs for a range of weak interactions, many applications in biochemistry benefit from the computational efficiency afforded by the GGA model. One practical approach to solve this problem and improve the GGA model at low computational cost involves the use of atom-centered dispersion-corrected potentials (von Lilienfeld et al., 2004). Benchmark studies show that the use of this empirical correction in combination with common DFT functionals can lead to large improvements in the performance of DFT for van der Waals forces (dispersion). In general, GGA functionals perform well for interactions involving electrostatics and, consequently, predict binding energies for ion–water interactions with relatively small errors (<2 kcal/mol) (Bucher et al., 2006). For hydrogen-bonded systems, the error in the energy of H bonds with common GGA functional is typically ~1 kcal/mol per H bond, or ~0.5 kcal/mol, with the atom-centered dispersion-corrected potential approach (Arey et al., 2009). Thus, in general, the fundamental interactions between water molecules or waters and ions can be described with reasonable accuracy at this level of theory.
To extend the applicability of AIMD to large biomolecules, one typically uses a QM/MM scheme. The basic idea of a QM/MM scheme is to describe part of the system at the quantum level, whereas the rest of the atoms (protein and solvent) is treated at the classical force field level. In this way, the computational resources can be concentrated on computing the electronic structure for a region of interest. Despite an important gain in computer efficiency with the QM/MM method, the length of the simulations remains a serious limitation of the method. QM/MM simulations can only be run for 10–100 ps on a system with roughly 500 QM atoms and 100,000 MM atoms. This renders the method impractical for applications that would require sampling on the nanosecond timescale. One important example is the calculation of converged free energies for the binding of an ion inside a channel. Classical studies performed more than 15 years ago (e.g., Roux and Karplus, 1993) showed that a large error (as large as ~5 kcal/mol) can occur in the computed binding free energy if sampling times are limited to ~50 ps. Thus, at the present time, computing converged free energies with AIMD remains a very challenging task and necessitates the use of enhanced sampling methods.
However, the field can be considered to be at its infancy, comparable to the state of classical MD 25 years ago. In the next decade, the method is likely to benefit strongly from the growth in computer power, the development of massively parallel computing, and the advent of linear scaling quantum mechanical methods. In particular for ion channel studies, it can be anticipated that this will allow for converged free energy calculations (on the nanosecond timescale) and lead to important new applications in ion channel modeling—among them, the accurate predictions of ion conduction and selectivity in ion channels from first principles.
Here, we will focus our attention on one current application of QM/MM simulations, which is the quantification of polarization effects. AIMD can be used as a tool to identify important challenges for force field developers. In particular, the method can measure the response of the electron clouds to moving ions (i.e., the ion-induced polarization), which is crucial to provide highly realistic models of conduction in ion channels. Polarization needs to be accounted for to increase the accuracy of classical MD simulations, but it often cannot be measured straightforwardly by experimental methods. Therefore, QM/MM simulations of ion conduction in ion channels can provide useful benchmarks that might assist the development of polarizable force fields.
In passing, we mention other reviews that have been written about the use of QM/MM simulations as a tool in the study of the enzymatic mechanisms (for reviews see Carloni et al, 2002; Rothlisberger and Carloni, 2006; Bucher et al., 2009a). The rate of enzymatic mechanisms can be measured experimentally and, compared with the prediction of the simulations, offer a true experimental validation of the method. Direct applications of QM/MM simulations to ion channel studies, on the other hand, have been very limited. Recent examples include the simulation of proton migration in aquaporin channels (Jensen et al., 2005) and the simulation of the coupling between ion conduction and the protonation states of pore-lining residues in the KcsA K channel (Bucher et al., 2007).
In the following paragraphs, we will discuss two common problems for nonpolarizable force fields that illustrate the need for better classical models: (1) the description of the coordination numbers of ions in K channels, and (2) the description of water polarization in the narrow pore of gramicidin (gA).
Classical MD simulations continue to play a central role in shaping our understanding of selectivity in ion channels. Luzhkov and Aqvist (2001) were the first to perform free energy perturbation calculations in the KcsA K channel to show that the narrowest region of the pore—the so-called selectivity filter—is selective for K+ versus Na+ in amounts comparable with experimental estimates. This finding was later confirmed by several groups and with a large variety of nonpolarizable force fields. We have also found that the intrinsic selectivity of the filter for K+ versus Na+ is still present when induced polarization effects are considered explicitly (Bucher et al., 2009b). Although one could argue that the selectivity of K channels may involve more than one region of the channel, and include multi-ion effects, the knowledge that the selectivity filter is intrinsically selective for K+ versus Na+ is interesting in itself. For example, it can come in useful for many practical applications in the field of nanotechnology, where such knowledge could be exploited to design artificial materials with tailored properties.
To a large extent, the selectivity of the channels can be understood by equilibrium thermodynamics, expressed through the relative free energies of different ions in aqueous solution versus the protein environment. However, some questions remain open. In particular, the nature of the physical interactions that are at the origin of the selectivity is still controversial. Initially, ion specificity of the selectivity filter has been explained as arising from the geometrical arrangement of ligands in the filter (Doyle et al., 1998), the so-called “snug-fit” hypothesis. However, MD simulations have shown that the selectivity filter is relatively flexible and that Na+ ions are able to interact favorably with the interior of the selectivity filter. A refined version of the snug-fit hypothesis was later proposed, in which the filter can adopt different conformations in the presence of K+ and Na+ ions (Lockless et al., 2007). In this model, the protein structure is expected to play an important role, and the property detected by the channel is the ion size. In addition, it has been proposed that the repulsion between carbonyl groups in the filter plays an important role in the selectivity (Noskov et al., 2004). The carbonyl dipoles have to come close to the Na+ ion to stabilize it inside the filter. However, this creates strong repulsions between the carbonyl dipoles. Finally, the selectivity has also been discussed as arising from the arrangement of the carbonyl groups in the filter, which can favor higher coordination for the ions compared with water (Bostick and Brooks, 2007; Thomas et al., 2007; Varma and Rempe, 2007). Central to this theory is the idea that the coordination number for each ion is not constant; it can vary in the protein environment.
Therefore, ion coordination numbers play an important role in the different theories of selectivity. However, because standard nonpolarizable force fields were mainly designed to describe average thermodynamic quantities, one may wonder if they can be used to reliably describe ion coordination numbers and fully address questions about the mechanisms of ion specificity. This seems especially doubtful in light of the fact that the close proximity of a positively charged ion is likely to strongly polarize the surrounding ligands, and ion-induced polarization effects can be expected to be particularly prominent for this case. Furthermore, an accurate determination of ionic coordination numbers is important because at the present time, information about the coordination numbers inside the channel is not directly available from experiments, and researchers have to rely solely on the results of MD simulations.
Recently, we used QM/MM simulations to describe the coordination numbers of K+ and Na+ ions in the KcsA channel (Bucher et al., 2010). In this approach, the selectivity filter was described with a first-principles (parameter-free) electronic structure method (DFT) that includes a description of induced polarization and charge transfer effects. Interestingly, the picture of the ion coordination given by QM/MM simulations differs somewhat from the one provided by classical MD based on nonpolarizable force fields. In particular, the ion coordination numbers inside the filter of the KcsA K channel predicted by the ab initio method are reduced by ~1. In addition, the coordination numbers for the two ions, in the filter and in water, are found to be very similar. We find that K+ displays an average coordination number between 6 and 7, both in water and in the filter. Similarly, Na+ prefers an average coordination number of ~5, both in water and in the filter. This result suggests that variation in the coordination numbers between the protein and the aqueous phase is small, and that a variation in the coordination numbers between the aqueous phase and the filter is not a requirement for selectivity. Furthermore, it supports the idea that the inclusion of polarization effects is indeed important to describe accurately the correct balance between ion–water and water–water interactions, which is crucial for discriminating between alternative theories of selectivity in ion channels.
The design of accurate microscopic models of ion channels is a long-sought goal of molecular modeling. Historically, the gA channel has been an important test case for microscopic models that aim to reproduce the current permeating through ion channels. It is a small bacterial channel, for which ample structural and functional information exist, allowing for an extensive comparison between experiments and simulations. The channel is formed by the assembly of two homo dimers, each made of 15 amino acids. It conducts monovalent cations and water molecules, which can permeate through the channel in a single file arrangement.
Rigid models using continuum electrostatic methods have been developed to model conduction in gA. However, in one study of the gA channel (Edwards et al., 2002), it was found that it is not possible to reproduce all the experimental data at this level of theory. Although one can fit the diffusion coefficient of ions so that the conductance of the channel is reproduced, in that case, the dependence of the conductance on the ionic concentration is not reproduced. Accurate MD simulations, on other hand, can potentially solve these problems. This realization provides an important motivation to develop further MD methods to accurately compute the conductance of ion channels.
However, at the present time, serious problems render the prediction of transport rates by MD simulations cumbersome. A first problem is the time span of a single ion permeation event (10–100 ns), which approaches the total sampling time available in a typical (force field–based) MD simulation. In most cases, enhanced sampling techniques (such as umbrella sampling or adaptive biased force methods) must be used to collect sufficient statistics to compute the free energy profile of ion conduction. In addition, a second problem is the accuracy of the underlying force fields, which in the case of nonpolarizable models, is limited by the mean-field description of polarization. It can be envisioned that in the next decade, faster computers will allow sampling times on the microsecond time scale, and that the most serious limitation of classical models will be the accuracy of the force fields.
In gA, the single ion potential of mean force for the translocation of a K+ ion has been computed with a variety of nonpolarizable force fields (Allen et al., 2006; Baştuğ, and Kuyucak, 2006). In these studies, the free energy profile of K+ conduction through the gA channel displays an energy barrier in the middle of the channel of ~10 kcal/mol that is inconsistent with the experimental observation that K+ conduction occurs at a rate close to the diffusion limit. The predicted current at low ion concentration is several orders of magnitude too slow because the current suppression is proportional to ~exp(−ΔG/kT), where ΔG is the barrier height for the single ion potential of mean force. A similar disagreement with experimental estimates has been observed for all standard nonpolarizable force fields (GROMOS, AMBER, and CHARMM), and it has been suggested that the lack of explicit polarization is at the origin of this discrepancy between experiments and classical MD (Allen et al., 2006; Baştuğ and Kuyucak, 2006).
Recently, we have used QM/MM simulations to study the impact of polarization on the conduction profile in the gA channel (Bucher and Kuyucak, 2009). This work was motivated by the need to quantify polarization effects on permeating waters because these effects are sometimes considered to be either negligibly small or, if sizeable, well described by a mean-field treatment. The QM/MM model of the gA channel was built so that the channel backbone atoms and the permeating waters and ions are described with a QM method, and the protein side chains, the membrane environment, and solvent are described at the force field level.
Special attention was given to the analysis of the dipole moments of water molecules in the hydration shell of K+. To compute water dipoles in different environments, we have used the formalism of maximally localized Wannier function. Wannier functions are a set of orbitals that are optimized using an iterative procedure to minimize their spatial spread, while still reproducing the correct electronic density for the system. The Wannier function are useful because they provide a description of the localization of electron pairs in orbitals that is intuitive for chemists, and include bonded orbitals and lone pairs. The center of these orbitals can be used to fit the charge parameters of a force field or to obtain rapid estimates of the impact of the environment on molecular dipoles in the condensed phase.
Not surprisingly, the polarization of waters was found to depend strongly on the local environment, in particular, on the presence of hydrogen bond acceptors and donors. When a K+ ion was placed at the center of the gA channel, high dipoles were obtained due to the induced polarization by the ion electrostatic field, and the single file alignment of water molecules. However, a different behavior of the water dipoles is observed in bulk, where cooperative effects between waters are absent, and the water dipoles are not enhanced by the presence of the K+ ion (Bucher and Kuyucak, 2008). This behavior of water dipoles is interesting because it may explain the discrepancy between experiments and the predictions of nonpolarizable force fields for the transport properties of K+ through the gA channel (Bucher and Kuyucak, 2009). Similarly, a recent study using the polarizable version of the CHARMM force field showed that good agreement between the computed free energy profile and experiments can be obtained in gA by including a description of polarization effects (Patel et al., 2009). Efforts are now being made to extract information from the QM/MM simulations of gA into polarizable force fields that can improve the description of ion conduction in narrow pores.
In our study of the bacterial KcsA K channel, we have attempted to extract optimized point charges for the residues of the selectivity filter using QM/MM potentials (Bucher et al., 2009b). We have found that the protein residues become strongly polarized in the presence of ions. In particular, the carbonyl groups in the selectivity filter of the KcsA channel interact directly with K+ ions, which leads to an enhancement of the carbonyl dipoles. The presence of sizeable ion-induced polarization effects suggests that these systems can benefit from the use of polarizable force fields. For example, different levels of ion-induced polarization for K+ and Na+ ions are likely to affect the description of K+/Na+ currents.
There is now growing awareness that ion channel studies can benefit from the development of polarizable force fields. Accurate microscopic models would be very useful, for example, in assisting the interpretation of electrophysiological data and to help understand the control of ion permeation by external factors, such as voltage, pH, concentration gradients, and drugs. Besides these factors, computational studies could, in principle, be used to predict and analyze the effects of various mutations on the ion channel conduction and ion specificity properties. This would be very useful to better understand genetic disorders that are caused by impaired ion channels, or to design better drugs that target ion channels.
An accurate description of polarization effects in ion channels will help close the gap between the predicted transport rates and experiments. In drug design, polarizable potential functions used to calculate drug target interactions are already producing some encouraging results, reaching chemical accuracy (<1 kcal/mol) in the calculations of binding free energies (e.g., see Khoruzhii et al., 2008).
Attention has now focused on improving current force fields by including polarization explicitly. To assist the development of polarizable force fields, it is important to single out systems of biological interest that could serve as defined performance benchmarks. Small ion channels offer useful benchmarks for force fields because the I-V curves can be determined experimentally at different ion concentrations and compared with computational predictions. By providing experimental data for many simple ion channels, the experimental community can play an important role in assisting the design of accurate force fields. In addition, in this Perspective, we propose to reinforce the importance of ab initio methods and in situ QM/MM data in parameterizing polarizable force fields.
In conclusion, we believe that QM/MM simulations will be highly useful in the future, both as a method to help parameterize polarizable force fields and as a method to study the conduction and selectivity properties of ion channels from first principles. At the present time, the short sampling times available (picosecond timescale) have been a serious limitation to the use of QM/MM simulations as a “production” tool in ion channel studies. However, in the next decade, such studies will begin to appear. Presently, the main use of the method in ion channel studies has been in assisting the parameterization of the next generation of force fields. In the condensed phase, QM/MM simulations should be used, when possible, to test and parameterize different functional forms. Historically, this represents an important paradigm shift because past force fields have been mainly required to describe average thermodynamics quantity and some structural and kinetic quantities. However, to incorporate explicitly electronic polarization into the force fields, it is clear that the parameterization process will need to rely more on microscopic quantities because experimental information about electronic polarization in biomolecular systems is scarce. Future collaborations between the experimental community and quantum chemists will be key to successfully design better force fields for ion channels.