Potassium channels are a family of membrane proteins that controls the rapid and selective conduction of K + and exclusion of Na + . 1 The exceptional properties of potassium channels and their importance in regulating cellular processes have stimulated decades of studies on their conduction mechanism. [1][2][3][4][5][6][7] The crystal structure of a model potassium channel, KcsA, 8 shows electron density in the four continuous ion binding sites in the selectivity lter (SF). The close distance of 3-3.5 Å between neighboring sites led to the hypothesis that K + only binds at two non-adjacent sites at a time to avoid strong electrostatic repulsion between ions. In this scheme, overlap between two alternating congurations, i.e. K + -water-K + -water and water-K +water-K + , produces the four peaks in the electron density. 8,9 This two-ion hypothesis, termed the "so knock-on" mechanism, has been used extensively in ion binding experiments and theoretical studies. Recently, however, this mechanism has been challenged by (i) a reanalysis of crystallography data showing the total occupancy is close to 4, 10 and (ii) long MD simulations of ion permeating events in which the rapid and selective K + conducting pathway contains no water. 10,11 Simulation studies to date have not fully claried the ion conduction mechanism. On one hand, the four-ion conguration was shown to be stable only in short simulations at 200 K, 10 likely the result of the slow dynamics at low temperature. On the other hand, two-ion and three-ion congurations accounted for more than 90% of long MD trajectories. 10,11 Also, it has been argued that simulation results are sensitive to force eld parameters, 12 as reected by other MD studies that support the so knock-on mechanism. [13][14][15] Furthermore, free energy simulations found that a direct knock-on mechanism and the so knock-on mechanism have comparable free energy barriers. 16,17 There have been numerous experimental studies supporting either mechanism. Crystallographic data are collected at low temperature, while other experiments are conducted at room temperature under either equilibrium or nonequilibrium conditions. Anomalous X-ray diffraction of NaK2K potassium channel showed high K + occupancy in the SF. 18 Solid-state NMR studies of NaK2K detected no water in the SF under physiological conditions. 19 2-D IR spectroscopy studies found that only the two-water two-ion congurations are compatible with the spectra of the SF, 20 while a later study suggested the 2D IR spectra could not differentiate the direct and so knock-on mechanisms. 11 Cuello and coworkers 21 showed that by removing the S3 site, a KcsA-G77A mutant stabilizes the water-K + -water-K + conguration. The G77A mutant has lower conductance but maintains K + /Na + selectivity, contrary to the previous conclusion that the direct knock-on mechanism is essential for selectivity. 11 The G77A mutant and the wildtype have almost identical structures except at S3, and similar ion binding affinities. Although the G77A mutant may not be directly relevant to wildtype KcsA, it is an interesting model system to study the ion conduction mechanism.
Due to the inconsistent results from various experiments and simulations, there is no consensus on the ion conduction mechanism of KcsA. The stability of ion congurations in the SF results from a competition among the interactions between ions, water and protein in the conned environment. 4,22 The polarization effect is signicant for such highly charged systems. [23][24][25] Polarization is the response of an electron cloud to the electric eld generated by its environment. The presence of ions dramatically changes the electric eld around the SF. So it is important to account for this effect by allowing protein and water to respond to this change in electric eld. Polarization can both enhance electrostatic interactions by increasing molecular dipole moments and reduce them through dielectric screening. For the interactions of ions, such as those in the SF, the screening effect dominates. [25][26][27] The screening effect can qualitatively affect ion binding. In fact, we showed previously that without the screening effect, Ca 2+ /Mg 2+ selectivity in most proteins will be inverted. [27][28][29] The importance of polarization for K + ion binding has also been recognized in recent studies. [30][31][32][33][34][35][36] Varma et al. highlighted the complexity of ion-protein interactions because polarization can both enhance and reduce the dipole moments of ion-coordinating molecules. 33 Rossi et al. found that the polarization of the threonine methyl groups at S4 site contribute to a few kcal mol À1 in ion binding energy. 34 Peng et al. and Ngo et al. showed that polarization signicantly reduces ion conduction barriers. 30,36 Lemkul and coworkers found that polarization is essential for the stability of G-quadruplexes. 37 G-quadruplexes are similar to the SF of potassium channels in that they contain a single line of K + ions with interatomic distances of $3.4 Å, and each binding site has 8 carbonyl groups surrounded by negatively charged groups. Despite these ndings, previous simulations of KcsA predominantly used xed-charge (nonpolarizable) force elds, 2 which could cause signicant errors for ion interactions. 26 Besides inaccuracy in force elds, the observation of both so and direct knock-on mechanisms in previous simulations could also be related to insufficient sampling. Transitions between congurations in the SF may require long simulations. Since the ion congurations and the structure of the SF are well dened, the relative stabilities can be conveniently calculated by free energy perturbation (FEP). 38 FEP utilizes alchemical pathways to calculate free energy differences between thermodynamic states or between molecules. FEP has been used widely for the study of protein-ligand binding, and is an efficient approach for small ligands with rigid binding pockets.
In this work, we endeavored to determine accurately the thermodynamic stabilities of various ion congurations through polarizable MD simulations and free energy calculations. Our results show that in contrast to previous simulations, the four-ion conguration in KcsA is thermodynamically stable and likely a frequent conformation during ion conduction. Additionally, we compared our simulations to the KcsA-G77A mutant data and experimental ion binding affinities.
Several hypotheses on ion occupancy have been proposed: full occupancy of 4 based on X-ray crystal structures, 10,18 reduced occupancy of 2-3 in the direct knock-on mechanism based on MD simulations, 10,11 and two-ion two-water conguration (Fig. 1A) based on alternative interpretation of X-ray structures and MD simulations. 9,21 Although the direct knock-on mechanism was considered to be consistent with full occupancy, MD simulations at room temperature have always produced reduced occupancy.
We predicted the relative stabilities of two-, three-and fourion congurations from binding free energies of K + /water and relative binding free energies between K + and water. Each absolute or relative binding free energy was calculated by FEP and MD simulations. The simulations mimic crystallographic conditions with a model system of KcsA solvated in water solution, while the effect of a lipid bilayer environment is also discussed. Two force elds that account for the polarization effect were chosen, the AMOEBA polarizable force eld and the CHARMM36m (C36m) force eld with electronic continuum correction (ECC). 25,26 For comparison, simulations with standard CHARMM were also included. AMOEBA explicitly represents polarization through induced dipoles on each atom. It has shown excellent accuracy for ion binding in both gas phase and condensed phase. 27,28,39 We further validated the accuracy of AMOEBA for ion-protein interactions against quantum mechanical (QM) methods (Fig. S1 andS2 †). ECC is a framework for modifying standard nonpolarizable force elds to represent polarization through a mean-eld approximation. In ECC, the charges of ionized groups and ions are scaled by sqrt(1/2) z 0.7, based on the argument that the electronic screening factor is about 2 for most organic media. 26 ECC has also been used in many MD simulation studies 25 thanks to its compatibility with popular MD soware and computational efficiency.
The results from both AMOEBA and C36m-ECC indicate that the four-ion conguration has the lowest standard free energy (Fig. 1B) despite the perceived strong repulsion between the ions. C36m-ECC simulations in DOPC lipid bilayer (Fig. 1B), and simulations with several variants of AMOEBA parameters (Fig. S3 †) also predict similar trends. The differences between the four-ion conguration and some three-ion congurations are moderate: under low K + concentration, reduced occupancy could become more stable. This result is consistent with the experimental observation of reduced ion occupancy in KcsA under low K + concentration of 3 mM. 8 In contrast, our simulations with the xed-charge CHARMM force eld favored the three-ion and two-ion congurations over the four ion-conguration, with 1,2,4-bound conguration being most stable. This result agrees with previous simulations that used xed-charge force elds. 10,11 Moreover, the general trend of the calculated free energies is insensitive to temperature (as low as 200 K). Even with the xed-charge force eld, the most stable, 1,2,4-bound conguration has a direct ion-ion contact at S1 and S2. This ion-ion contact occurs because the four negatively charged D80 residues near S1, in addition to the strong dipole moments of the carbonyl groups, favor ion binding at S1 and S2 and the articially strong ion-ion repulsion prevents binding of three continuous ions. From the comparison between various force elds, especially between C36m-ECC and C36m where the only difference is effective polarization of ions, it can be seen that polarization strongly inuences ion binding thermodynamics. Only simulations with polarization are consistent with the full ion occupancy observed in crystal structures. Despite similar trends regarding the stability of the four-ion conguration, some noticeable differences between AMOEBA and C36m-ECC exist (Fig. 1B). First, the water-K + -K + -K + conguration is most stable among 3-ion congurations in AMOEBA, while several 3-ion congurations have similar free energies as predicted by C36m-ECC. Second, the free energy difference between water-K + -water-K + and water-K + -vacancy-K + , i.e. the water binding free energy for this conguration, varies signicantly. In both AMOEBA and C36m results, the water binding free energy is close to 0, meaning the congurations with and without water at S3 are both possible. With C36m-ECC, the water binding free energy is about À4 kcal mol À1 , which seems unusual.
The stability of the four-ion conguration is largely enthalpydriven (Fig. 1C). The polarization contribution to enthalpy favors the four-ion conguration. The variations in polarization, electrostatics and vdW enthalpies are much larger than the variation in the total enthalpy, indicating compensation between these energy components.
The free energy results were conrmed by long MD simulations of KcsA in lipid bilayer. With effective polarization, starting from the four-ion conguration, the rst ion leaves the SF and the h ion enters the SF, resulting in a four-ion conguration stable in 500 ns simulations. Without polarization, the SF quickly adopts 1,2,4-bound conguration and remains stable through 500 ns simulations. This result is consistent with our free energy calculation. However, the nal structures also depend on the initial structure, which indicates insufficient sampling. Starting from the water-K + -water-K + conguration, long-lived congurations K + -K + -water-K + and K + -water-K +vacancy-K + (the rst K + at the entry of S1) were observed from simulation with and without effective polarization, respectively.
Multiple binding/unbinding of K + at S1 were observed in two 100 ns AMOEBA simulations of KcsA in 0.15 mM KCl solution. The four-ion conguration accounts for 61% of the last 80 ns simulations, with the rest being three-ion congurations (Fig. S5 †). This percentage is converted to a standard binding free energy for the fourth ion of À1.4 AE 0.9 kcal mol À1 , in excellent agreement with À1.3 AE 0.6 kcal mol À1 from the free energy calculation (Table 1 and Fig. 1B). AMOEBA simulations of KcsA in DOPC bilayer were stuck in the water-K + -K + -K + conguration (Fig. S5 †), which is also predicted to be a low free energy conguration in water solution (Fig. 1B).
To better understand the polarization effect on ion binding and ion-ion repulsion, the energy components as a function of ion-ion distance were analyzed. Starting from the crystal structure of KcsA, the rst K + ion was gradually moved out of the SF (Fig. 2). AMOEBA and C36m-ECC both predict that structures with the rst K + inside and outside the SF have comparable energies separated by a barrier, while C36m and AMOEBA without polarization predict the rst K + will be repelled out of the SF with almost no barrier (Fig. 2A). Notably, the barrierless repulsion is an integral part of the previous direct knock-on mechanism. 10 At short ion-ion distance between 3 and 5 Å, the electrostatic interaction is repulsive, while the polarization energy is nearly proportional to the opposite of electrostatic energy, therefore signicantly reducing ion-ion repulsion. The polarization comes from the surroundings of the ions. It is long ranged, not fully converging aer 10 Å (Fig. 2B). Previous work also pointed out that models that explicitly include the protein environment are more accurate than active-site models. 40 The reason why polarization effectively reduces ion-ion repulsion is as follows: when the two ions come closer, the combined electric eld produced by the two ions also becomes stronger, leading to a much more favorable polarization energy which is proportional to the square of the electric eld. This result is a classic example of dielectric screening. If the environment is treated as an electronic continuum, the screening effect can be modeled by a dielectric constant or scaled effective charges. 26,[41][42][43][44] In addition, because this screening effect mostly acts on ion-ion repulsion, it cannot be straightforwardly modeled by modication of ion-protein vdW interactions as in previous work. 45 Ion conduction is driven by a membrane potential and/or ion concentration gradient, therefore it is important to consider the effect of membrane potential on the conduction mechanism. The membrane potential of neurons lies between À70 and 40 mV. The holding potential used in electrophysiological experiments is up to 200 mV. 1 The effect of the membrane potential was studied by using the model system in Fig. 2 with either an external electric eld or extra ions on each side of the membrane. As shown in Fig. S6, † a membrane potential of 500 mV only shis the relative energy curve in Fig. 2A by 0.8 kcal mol À1 , which is much smaller than the barrier and the effect of polarization. Therefore, the membrane potential has a negligible effect on short-range ion-ion interactions. However, the membrane potential is important for MD simulations of ion permeation since it directly affects ion conduction rate and, in some cases, the conformation of ion channels.
One recent piece of evidence for the so knock-on mechanism comes from the observation of 2,4-bound conguration in KcsA-G77A (Fig. 3A) and KcsA-G77C mutants. 21 The mutants have abolished binding at S3 due to the rotation of the backbones of V76. The G77A mutant is a K + -selective channel with a reduced conductance $32 times lower than that of the KcsA wildtype (KcsA-WT). The G77A mutant also has a similar apparent K + binding affinity as KcsA-WT. 21 Although the mutants could have a different conduction mechanism than the KcsA-WT, they serve as model systems to study how conduction mechanism is affected by the structure of the SF.
It was argued that if KcsA-WT has four ions in the SF, KcsA-G77A should have three ions at S1, S2 and S4 (Fig. 3A). 21 As polarizable force elds produce weak ion-ion repulsion in the SF, it is interesting to see whether the 2,4-bound conguration in KcsA-G77A can be reproduced. The free energy calculation of KcsA-G77A shows that indeed the 2,4-bound K + conguration is more stable than the 1,2,4-bound conguration and the one-ion conguration (Fig. 3B). Although the mutation only affects S2 and S3, both S1 and S3 are free of K + ions. The coupling between S1 and S3 is due to the interaction between ions: K + at S2 is close to an in-plane binding position between S1 and S2 (Fig. 3A), which excludes K + at S1 despite the reduced ion-ion repulsion predicted by polarizable force elds. Previously, the coupling between ion binding at S1 and S3 or between S2 and S4 has been considered as evidence of the two alternating congurations in the so knock-on mechanism. 21,46 Apparently, the coupling of reduced occupancy at S1 and S3 in KcsA-G77A does not contradict the full ion occupancy in KcsA-WT.
Ion binding experiments provide high-quality thermodynamics data that can be used to validate force elds. Ion binding in potassium channels involves multiple binding events, competition between ions and conformational changes. Thus, it has been difficult to interpret the experimental ion binding data. Lockless and coworkers 47 showed that K + binding in KcsA and MthK potassium channels in isothermal titration calorimetry (ITC) conditions is the competitive binding between Na + and K + . Using ITC with a competitive binding model and several KcsA mutants that only adopt one SF conformation, the intrinsic binding affinities for different conformations were determined. The K + binding affinity was found to be $0.1 mM, and the Na + / K + selectivity varies between 10 2 to 10 4 for different channels and conformations. A much stronger K + binding affinity of $2 mM was determined by thermal denaturation experiments. 48 For the collapsed conformation of KcsA, since the SF is unstable The two values from top to bottom are free energy changes from 2 K + to 3 K + and 3 K + to 4 K + in the SF, respectively. 3 K + free energy is calculated by sum of partition functions of all four congurations. 2 K + free energy is from the water-K + -water-K + conguration. b Free energy change from 2 Na + to 1 Na + /1 K + in the SF. c Free energy change from 1 Na + /3 K + to 4 K + in the SF. d Experimental binding free energy was calculated from DG ¼ ÀRT ln [K d /(mol L À1 )], with K d (K + ) ¼ 67 and 350 mM, K d (Na + ) ¼ 9 and 0.48 mM in wildtype (conductive) and M96V mutant (collapsed) KcsA, respectively. 47 without ions 48 and there are a maximum of two K + , the binding affinity likely corresponds to binding of a second K + . The conductive conformation of KcsA is observed at higher K + concentration than the collapsed conformation, so it should contain at least 2 K + , and the binding affinity may be related to the binding of the third or the fourth K + . The binding free energy for each conformation can be calculated separately by simulations since the simulation time scale is relatively short.
By integrating the free energies of various congurations (eqn (S1), Table S1 †), we calculated the total free energy for one-ion to four-ion congurations and derived the binding free energy for each binding event (Table 1). The calculated binding free energy for a third K + in the KcsA conductive conformation is À6.8 kcal mol À1 , close to the ITC value of À5.6 kcal mol À1 , 47 while the binding of a fourth K + seems too weak (À1.3 kcal mol À1 ) to be measurable. The calculated K + binding free energy in the collapsed conformation of KcsA (À4.9 kcal mol À1 ) and Na + /K + relative binding free energies (À1.6 to À2.0 kcal mol À1 ) also agree well with experiments 47 (À4.7 and À2.1 to À2.9 kcal mol À1 , respectively). One intriguing question from the KcsA-G77A mutant experiments is why KcsA-G77A and KcsA-WT have similar binding affinities (0.30 and 0.29 mM) 21 if the ion congurations are different. KcsA-WT adopts the collapsed conformation at low K + concentration and transitions to the conductive conformation when K + concentration is above 5 mM, although it is unclear whether the binding takes place before or aer the conformational change. 47,49 KcsA-G77A has the conductive conformation at both low and high K + concentrations. In our free energy calculations, the binding free energy for KcsA-G77A (À6.0 kcal mol À1 ) is close to both the collapsed (À4.9 kcal mol À1 ) and conductive conformations (À6.8 kcal mol À1 ) of KcsA-WT, considering the accuracy of MD The polarization energy for K + and its environment within a certain distance is also computed. "Pol (K + )" means polarization energy for four K + in the SF, "Pol (5 Å)" means the polarization energy for K + and atoms that are within 5 Å of K + at any point during the distance scan, and likewise for other polarization energy. The energy was calculated using crystal structure (PDB ID: 1k4c) embedded in DOPC bilayer, except that the coordinate of K + at S1 was modified and two K + ions above S1 were removed.
Fig. 3 Relative free energies of ion configurations in KcsA-G77A mutant. (A) Crystal structure of KcsA-G77A (PDB: 6nfu) and a hypothetical 1,2,4-bound configuration from MD simulations; (B) relative free energies of ion configurations calculated from MD simulations. Blue and green circles in the cartoon representations stand for K + and water, respectively. The uncertainties of the free energies are smaller than the size of the symbols.
simulations (Table 1). This data suggests that the binding affinities could be similar despite different conformations or ion congurations.
Since the SF is saturated by K + at equilibrium, the most likely conduction mechanism is the single vacancy mechanism, 50 where one vacancy is created by unbinding of the rst/last ion and then the vacancy is gradually transferred along the SF. The so knock-on mechanism, which proceeds though alternating congurations of water-K + -water-K + and K + -water-K + -water, 9,13 is unfavorable due to the high free energy of the two-ion conguration (Fig. 1B), unless it has much lower barrier than the single-vacancy mechanism. To assess the barriers of the single vacancy and the so knock-on mechanisms, we calculated the potential of mean force (PMF) as a function of ion coordinates by umbrella sampling. 51,52 The intermediate steps in the single vacancy mechanism and the so knock-on mechanism are illustrated in Fig. 4A andC, respectively. Previous simulations of the so knock-on mechanism found K +water-K + -K + as an intermediate state. 13 This structure was also sampled in our simulations; it is slightly unfavorable compared to the water-K + -vacancy-water-K + conguration (indicated by (iii) in Fig. 4B andC), but does not affect the overall barrier. According to the PMF (Fig. 4A), it is unfavorable for the four-ion conguration to lose one ion at S1 or S4, consistent with the free energy results in Fig. 1. The overall barrier for the single vacancy mechanism is $9.5 kcal mol À1 and the highest barrier for individual steps is $6.5 kcal mol À1 (Fig. 4A). The barrier is considerably larger than 2-3 kcal mol À1 reported in previous simulations. 13,16 The large barrier is due to the close carbonyl O-O distances in simulations, which in turn is caused by the backbone torsional angles, ion-protein interactions, and the restraints from the protein scaffold. The distances between in-plane diagonal carbonyl oxygen atoms of T75 and V76 in the crystal structures are 4.49 and 4.60 Å. The corresponding distances in equilibrium MD simulations are 4.31 and 4.46 Å. There is also a noticeable difference in the V76 f angle (À78.6 in MD and À68.1 in the crystal structure). The equilibrium K + -O distance in K + -dialanine dimer is 2.55 Å (Fig. S2 †). So the O-O distance needs to increase to $5 Å during ion conduction. A stiff protein scaffold or backbone torsions will lead to a large barrier. Fine tuning of the force eld parameters is necessary for realistic simulations of the ion conduction. The equilibrium ion binding properties are less affected because the K + -O distance at the binding site is larger than at the in-plane position. Nevertheless, the barrier of the single vacancy mechanism is lower than the barrier of the so knock-on mechanism ($12.0 kcal mol À1 , Fig. 4B andC). The higher barrier for the so knock-on mechanism can be explained by the stronger ion binding in the half-lled conguration and thus the larger energy cost for ions leaving their favorable binding pose.
The difficulty in inferring ion conduction mechanism of potassium channels mainly comes from the contradiction between strong ion-ion repulsion and the observed ion occupancy at four continuous binding sites in the SF. Initially, the so knock-on mechanism was proposed to resolve this issue by hypothesizing that ions are separated by one water molecule. 9 The so knock-on mechanism requires two energetically balanced congurations, K + -water-K + -water and water-K +water-K + . In the direct knock-on mechanism, ions transiently occupy three continuous sites and the repulsion from the rst two ions helps the third ion overcome the barrier. 10,11 The direct knock-on mechanism explains the Na + /K + selectivity. However, it contains at most three ions in the SF simultaneously. To explain the ion occupancy in crystal structures, it was argued that the four-ion conguration only exists at very low temperature, which was not veried by converged simulations. Therefore, neither mechanism gives a satisfactory explanation of the crystal structure.
Although ions of like charge strongly repel each other, the SF also has a large negative potential that could possibly compensate for the ion-ion repulsion. To determine the balance between the two forces, accurate electrostatic models are required. Polarization effects are critical for highly charged systems, but they were oen neglected in previous simulations of ion channels. Our data show that polarization effectively reduces the ion-ion repulsion through dielectric screening. As a result, the four-ion conguration observed in crystal structures is thermodynamically stable. Considering the thermodynamic stabilities, the most likely conduction mechanism is through reversible binding of ions and the movement of one ion vacancy, or the single-vacancy mechanism. This mechanism is much simplied and consistent with crystal structures.
The single vacancy mechanism suggested by our data is close to the direct knock-on mechanism from previous work, but there are some noticeable differences. In the single vacancy mechanism, the four-ion conguration is stable, and the onevacancy states incur a small free energy cost. In the previous direct knock-on mechanism, three-ion states are most stable, and the conduction also involves four-ion and two-ion states (see, e.g., Movie (S1) in ref. 10). This difference could affect the conduction rate and the conformation of the selectivity lter, which can be veried by comparing to electrophysiological 50 and spectroscopic experiments. 20 Additionally, previous simulation results were sensitive to force eld parameters; slight changes in the parameters could lead to the so knock-on mechanism. 12 In contrast, our prediction of full ion occupancy is insensitive to force eld parameters once polarization is included. Our hypothesis is mainly based on equilibrium ion binding properties using the crystal structures, while actual ion conduction conformations may be different. Future simulations under realistic conditions and at larger time scale are needed to establish the detailed conduction mechanism.
This work has clearly explained the ion occupancy in crystal structures and ITC binding affinities, both of which are directly related to ion binding thermodynamics. In addition, ITC experiments of the MthK K + channel found that there are two K + binding events, each with a Hill coefficient of $2 for Na + displacement, which is only possible if there are at least four Na + ions before K + binding. 53 Recently, the KcsA-G77A mutant was found to have a stabilized 2,4-bound ion conguration characteristic of the so knock-on mechanism. 21 Our data shows that this is because it has a different ion binding property from the wildtype. There have been other experimental studies contradicting the direct knock-on mechanism, and they have been explained partly in recent work of de Groot and coworkers. 10,11 The ion-water coupling ratio derived from streaming potential measurement of KcsA was 1.0, 54 but the applied osmotic pressure could lead to an ion-depleted, waterpermeable state, which is distinct from the conductive state. 10 2D IR spectroscopy combined with MD simulations found that only a combination of water-separated two-ion congurations, including a structure with one ipped carbonyl group, was compatible with the spectra, 20 while it was shown later that the spectra could be reproduced by using only congurations with direct ion-ion contact and the carbonyl-ipped states were unnecessary. 11 As most of the structures with direct ion-ion contact also exist in the single vacancy mechanism, the 2D IR study is also compatible with our data.
MD simulations have been valuable for the study of both ion conduction and activation mechanisms of potassium channels. [55][56][57][58] Considering the importance of polarization for ion binding thermodynamics, we strongly advocate the use of polarizable force elds (such as AMOEBA, CHARMM Drude 59 ) or force elds with effective polarization for ion channel simulations. The charge scaling method (ECC) is a viable approach to effective polarization when computational cost and/or soware compatibility are a concern. We have shown that a scaling factor of 0.7 yields the same ion occupancy as explicit polarizable force elds, but deviation exists for the relative stabilities between various congurations. Fine tuning of the scaling factor and possibly vdW parameters 25 may be needed for more realistic simulations. Alternative approaches are QM-based methods that properly account for both local and long-range interactions, such as QM/MM 40 and quasi-chemical theory. 6,35
Systems with KcsA WT in collapsed, partially open 60 and closed conductive conformations and G77A mutant (PDB codes: 1k4c, 1k4d, 36 and 6nfu) in water boxes or embedded in DOPC bilayers (Fig. S4 †) were set up by using the CHARMM-GUI. 61 The box sizes for the water and DOPC systems were about (72 Å) 3 and 76 Â 76 Â 92 Å, 3 respectively. Salts were added to give a concentration of 0.15 M NaCl for the water box to mimic the crystallography conditions and 0.4 M KCl for the DOPC bilayer systems to mimic the high local concentration during ion conduction. For the DOPC system, a multistep procedure recommended by CHARMM-GUI was used to for initial relaxation.
The AMOEBA force eld parameters were validated by comparing to QM and experimental ion binding energies (Fig. S1 andS2 †). The accuracy of our force eld is comparable to popular density functional theory methods.
Simulations with the AMOEBA force eld were performed using Tinker-OpenMM. 62,63 The AMOEBA 2013 parameters for protein, 64 the ion parameters developed by Wang 65 and the DOPC parameters developed by Li and coworkers 66 were used. Larger polarizabilities for carbonyl groups from Liu et al. 67 were used since they better describe ionic interactions (see ESI †). Simulations with CHARMM36m force eld 68 were performed using GROMACS 2018. 69 The water box and DOPC box were used for binding free energy calculations and direct MD simulations, and the DOPC box was used for PMF calculation with AMOEBA and direct MD simulations. The Glu71 sidechains were protonated according to previous studies, 11 and all other residues were assigned default protonation states. The collapsed and closed conductive conformations are used for binding free energy calculations and the partially open conformation was used for PMF calculation. To keep the proteins at their starting conformations, 23 at-bottom position restraints were applied to all alpha carbons except those in the SF (residues 74 to 81). The restraints have a force constant of 5 kcal mol À1 ÅÀ1 for deviations larger than 2 Å.
Long MD simulations without ion restraints were performed, including two 500 ns simulations with DOPC box for each combination of CHARMM/CHARMM-ECC and four-ion/two-ion starting structure, and two 100 ns AMOEBA simulations with either water box (0.15 mM KCl) or DOPC box staring from the four-ion conguration.
Free energy perturbation (FEP) and the double decoupling method 70,71 was used to calculate the standard binding free energies in various congurations to derive the relative free energy between different congurations. For example, the DG between two congurations K + -vacancy-K + -K + and K + -K + -K + -K + was calculated as DG for moving K + from gas phase to S2 minus DG for moving K + from gas phase to the solution. The relative free energy between Na + -vacancy-vacancy-Na + and Na + -vacancy-vacancy-K + was calculated as DG for mutating Na + at S4 to K + minus DG for mutating an aqueous Na + to K + . The list of FEP calculations is tabulated in Tables S2 andS4. † Each free energy change was calculated through a series of alchemical transitions. For systems with net charge, there is a nite-size effect that depends on not only the box size but also the distribution of charges in the system. To remove such nite-size effect from nal DG values, the same or similar systems (which differ by four G77A mutations) were used for DG in the aqueous phase and DG in the SF. Alternatively, the nite size effect can be corrected analytically or by keeping the system charge-neutral during the alchemical transition. Full details of the calculation can be found in the ESI. † In the ECC approach, the solvation free energy consists of a nuclear contribution and an electronic contribution. 72 The nuclear contribution is calculated by standard FEP using the scaled charge. The electronic contribution can be calculated by the effective Born radius and the electronic dielectric constant, which does not depend on the atomic congurations. Therefore, the electronic contribution is canceled out in the nal binding free energy and is omitted in this work. It was shown that traditional nonpolarizable force elds perform well for high-dielectric environment such as aqueous solution because the missing electronic contribution is almost exactly captured by the overestimation of the nuclear contribution, but they perform poorly for low-dielectric environment. 72
The crystal structure of KcsA (PDB code: 1k4c) was used as a model structure to understand the effect of polarization on ion-ion interaction. The initial structure prepared for free energy calculation of KcsA embedded in DOPC bilayer was rst cleaned by energy minimization. Then the two K + ions above S1 were removed. The distance between K + ions at S1 and S2 was varied by modifying the z-coordinate of K + at S1. The vdW cutoff and Ewald parameters were same as those in free energy calculation. The contribution of each atom to the polarization energy was calculated by the ANALYZE program in Tinker 7. Calculation with an external electric eld was performed using a modied version of Tinker 8. 63 Ion permeation potential of mean force (PMF) calculation
The ion conduction barriers as a function of the z-coordinates of the ions were calculated by 1-D or 2-D WHAM methods using AMOEBA. For the direct knock-on/single vacancy model, each movement of the vacancy involves a different ion, and requires a separate 1-D WHAM calculation. For the so-knock on mechanism, the positions of the two ions initially at S4 and the cavity were used as reaction coordinates for 2-D WHAM calculation. WHAM consists of simulations in overlapping windows, where reaction coordinates (ion positions in this work) were restrained at different values in each window. This ensures even sampling of the reaction coordinates and helps overcome high barriers. The probability distributions of reaction coordinates in each window are then combined to reconstruct the free energy prole along the reaction coordinates. The WHAM code 73 was used for the analysis. Materials and Methods Quantum mechanical (QM) calculations. QM gas-phase binding and interaction energies were calculated to validate the AMOEBA force field. Geometry optimization was carried out with MP2/aug-cc-pVTZ/def2TZVPP by using Gaussian09 1 or Psi4 package. 2 MP2, DFT and CCSD(T) single-point calculations were performed using Psi4. The QM energies were compared with AMOEBA with several sets of parameters, [3][4][5] as explained in Fig. S1 caption. Counterpoise correction in Gaussian09 and Psi4 was applied for all interaction energies. MD simulations. For AMOEBA simulations, the equation of motion was integrated by the RESPA integrator with an outer timestep of 2 fs and the temperature was controlled by the Bussi thermostat at 298 K. The electrostatics was treated by particle-mesh Ewald (PME) method with a real-space cutoff of 8 Å and grid size of 0.9 Å, and the polarization was solved by the OPT4 algorithm. The vdW cutoff was 12 Å. For CHARMM simulations, the leapfrog integrator was used and all hydrogen atoms were constrained by SETTLE or LINCS algorithms to allow for a timestep of 2 fs. The Berendsen/Bussi thermostat for equilibration/production and the Berendsen barostat for equilibration were used to maintain the system temperature and pressure. The cutoffs for realspace electrostatics and vdW were set to 12 Å and long-range electrostatics was treated by PME with grid size of 1.2 Å. The initial equilibration of the DOPC systems employed semi-isotropic barostat. The DOPC systems were equilibrated through the procedure recommended by CHARMM-GUI, which consists of multiple steps with decreasing restraints on protein backbone and sidechain heavy atoms and on the torsional angles of the lipids.
Free energy calculation. The standard protocol as in our previous work [6][7] was used to calculate the standard binding free energies. In absolute binding free energy calculations using AMOEBA, a total of 18 alchemical states were set up to connect the two end states. The electrostatic interactions between the ligand (ion or water in our simulations) and the environment were first gradually decoupled in 10 steps (λ = 1.0, 0.9, 0.8, …, 0.0) before the vdW interactions were turned off in 8 steps (λ = 1.0, 0.9, 0.8, 0.7, 0.6, 0.55, 0.5, 0.0). The default vdW soft-core parameters in Tinker were used. A restraint on the ligand in the decoupled state was applied to avoid bad convergence and an analytical correction was added to account for the free energy change between the standard state (1 mol/L in both gas phase and solution phase) and the constrained state. Flat-bottom restraints between the ligand and the center of mass of the binding sites defined by carbonyl groups were used to maintain the designated ion configurations in the bound state. These flat-bottom restraints have a radius of 1.8 Å and a force constant of 50 kcal/mol/Å 2 , which specify a volume similar to the size of the binding sites. The restraints were gradually changed to a harmonic restraint with force constant 15 kcal/mol/ Å in the fully decoupled state (λ = 0.0). The free energy change from the decoupled state to the standard state in gas phase was analytically calculated to be 6.23 kcal/mol. 8 For each state, the simulations consist of 1-ns NVT equilibration and 2-ns NVT production. The production simulation was used for free energy perturbation using the Bennet acceptance ratio (BAR) method.
A similar protocol was used for GROMACS CHARMM calculations, with a few differences noted below. A total of 20 alchemical states were set up, with 10 states for electrostatics and 10 states for vdW (λ = 1.0, 0.9, 0.8, 0.7, 0.6, 0.55, 0.5, 0.4, 0.2, 0.0), due to the different vdW soft-core function. The flat-bottom restraints were kept constant with a radius of 1.8 Å and a force constant of 47.8 kcal/mol/ Å 2 (or 40000 kJ/mol/nm 2 in GROMACS units and convention), because mutation of the radius was not supported. The analytical standard state correction was 2.49 kcal/mol. The simulations consist of 2-ns NPT equilibration and 4-ns NVT production. VdW softcore parameters were sc-alpha = 0.5, sc-power = 1, sc-r-power = 6, sc-sigma = 0.3. No coulomb softcore was used.
The AMOEBA relative binding free energies were calculated by mutating the force field parameters. 15 steps were used for water-K+ relative binding free energy. The O vdW/polarizability parameters were mutated to those of K+ in 3 steps, then the O charge was changed to +1 in 10 steps while O dipole/quadruple and H multipole/polarizability were changed to 0. Last, the H vdW was turned off in 2 steps. For Na-K relative binding free energy, all parameters were linearly changed from Na to K in 4 steps. 1-ns equilibration and 2-ns production simulations were conducted for each state.
Using the 4-ion configuration as a reference, the free energy of each 3-ion configuration was derived from one double decoupling calculation. In addition, the free energy of water-K+-K+-K+ was also calculated by water-K+ relative binding at S1 using AMOEBA, which gives 0.65 ± 0.16 kcal/mol compared to 0.76 ± 0.21 kcal/mol from K+ binding at S1. For AMOEBA, the free energy of water-K+-water-K+ was calculated by water-K+ relative binding at S3. Then water-K+-vacancy-K+ was calculated by additional water binding at S3. The free energy of water-K+-vacancy-K+ from this path is 5.06 ± 0.30 kcal/mol, compared to 5.60 ± 0.30 kcal/mol from K+ binding at S3. The agreement from different paths further verifies the convergence of the free energy calculation. For CHARMM, water-K+-vacancy-K+ was first calculated by K+ binding at S3, and then water-K+water-K+ was calculated by an additional water binding at S3.
For the collapsed conformation of KcsA and the conductive conformation of KcsA-G77A, the 2ion configurations in the crystal structures were used a reference.
The AMOEBA absolute binding free energy calculations used the amoebapro13 parameters. The results of mod1 and mo2 parameters were obtained by FEP from amoebapro13 to the modified parameters. Since both are small modifications, three steps were enough to obtain converged free energies.
The total free energy for 1-ion, 2-ion or 3-ion configurations can be calculated by summing up the partition function of each configuration,
(S1)
When one configuration has much lower free energy than others, the total free energy can be approximated by the free energy of this configuration. The two 2-ion configurations that we calculated are assumed to have the lowest free energies. Using Eq. (S1) and the data from Table S1, the binding free energy for the second, third and fourth ion can be calculated and compared to experiment. For G77A mutant, we were unable to calculate the free energy for the only one ion at S4 since this configuration is very unstable. Similarly, the Na+-K+ relative binding free energy can be calculated by integrating possible configurations. For the conductive conformation of KcsA, only Na+-K+-K+-K+ and K+-K+-K+-Na+ were calculated and they have similar free energy. For the collapsed conformation of KcsA and the conductive conformation of G77A, the relative free energy was calculated as ΔG from 2 Na+ configuration to 1 Na+/1 K+ configurations. The most stable 1 Na+/1 K+ configuration are Na+vacancy-vacancy-K+ and water-Na+-water-K+, respectively.
The relative enthalpy of each ion configuration was calculated by the average energy from 2 ns NVT simulations of the same system at respective configuration. The same restraints as in the FEP calculation were used to maintain the ion configurations.
In both 1-D and 2-D PMF calculations, the spacing between windows was 0.35 Å and the force constant was 20 kcal/mol/Å. The initial structures were generated by pulling the ion(s) to the desired position(s) in 200-500 ps simulations, and at least two 200-ps simulations were performed at each window in the production run. The center of mass of the SF backbone atoms were restrained at the starting position. The 1-D PMF of the single vacancy model consists of multiple steps. First the ion at S1 was pulled out of the SF to the extracellular side, then the ion at S2 was pulled to S1, then the ion at S3 was pulled to S2, then the ion at S4 was pulled to S3, and in the last step the ion in the cavity was pulled into S4. For the 2-D PMF of the soft knock-on mechanism, the z-coordinates of the ions at S4 and the cavity were used as the reaction coordinates. The ion at S4 was moved upward to S2, and the z-direction distance between the two ions were scanned between 3.0 and 7.5 Å, which covers most likely intermediate states.
Comparison of experimental, QM and AMOEBA gas-phase binding energy. The binding energy is evaluated at 0 K. C: threshold collision-induced dissociation; E: single temperature equilibrium; K: kinetic method results. 9 The geometries and deformation energies are obtained by using MP2/aug-cc-pvtz. CCSD(T) interaction energies are calculated by MP2/CBS + δCCSD(T)/aug-cc-pvdz. The aug-cc-pvdz basis set was used with B3LYP-D3. The parameter "mod1" is amoebapro13 + modification of C=O dipole based on ion-dipeptide interaction; the parameter "mod2" is amoebapro13 + modification of C=O polarizabilities according to Liu et al. 3 The parameter "mod2" is used in simulations reported in the main text. a In the "Conf." column, the four-letter code indicates the species at the four binding site. "K", "w" and "0" indicates K+, water and vacancy, respectively. For example, "wK0K" means water-K+-vacancy-K+. Also, vacancy at the 1 st or 4 th position means water-bound and is only for convenience, since no restraints were applied to prevent water from occupying S1 and S4. b "amoeba13 → mod1" means the relative free energy for changing the force field parameter from amoeba13 to mod1, and similarly for "amoeba13 → mod2". c "sd" means one standard deviation for the value in the preceding column.
© 2021 The Author(s). Published by the Royal Society of Chemistry
© 2021 The Author(s). Published by the Royal Society of Chemistry Chem. Sci., 2021, 12, 8920-8930 | 8925
© The Author(s). Published by the Royal Society of Chemistry Chem. Sci., 2021, 12, 8920-8930 | 8929