The selectivity of the Na+/K+-pump is controlled by binding site protonation and self-correcting occlusion

The Na+/K+-pump maintains the physiological K+ and Na+ electrochemical gradients across the cell membrane. It operates via an 'alternating-access' mechanism, making iterative transitions between inward-facing (E1) and outward-facing (E2) conformations. Although the general features of the transport cycle are known, the detailed physicochemical factors governing the binding site selectivity remain mysterious. Free energy molecular dynamics simulations show that the ion binding sites switch their binding specificity in E1 and E2. This is accompanied by small structural arrangements and changes in protonation states of the coordinating residues. Additional computations on structural models of the intermediate states along the conformational transition pathway reveal that the free energy barrier toward the occlusion step is considerably increased when the wrong type of ion is loaded into the binding pocket, prohibiting the pump cycle from proceeding forward. This self-correcting mechanism strengthens the overall transport selectivity and protects the stoichiometry of the pump cycle. DOI: http://dx.doi.org/10.7554/eLife.16616.001


Introduction
The Na + /K + -pump is a primary active membrane transporter present in nearly all animal cells. It belongs to the P-type ATPase family, which utilizes the energy released from ATP hydrolysis to move ions against their concentration gradients across a membrane barrier. The ion species transported by the pump under physiological conditions are Na + and K + . Like many other membrane transporters, the Na + /K + -pump works according to an 'alternating-access' ion transport mechanism, with the bound ions accessible from only one side of the membrane at a time. The consensus scheme of the pump cycle is known as the 'Albers-Post' cycle (Albers, 1967;Post et al., 1969). It involves two major conformations, E 1 and E 2 , with inward-and outward-facing ion binding sites, respectively. In each cycle, the E 1 conformation binds three Na + from the cytosol and exports them using the energy from ATP hydrolysis. After external release of Na + , binding of extracellular K + promotes the structural transition to the occluded E 2 state, which finally imports two K + as binding of ATP returns the pump to the E 1 conformation. E 1 /E 2 conformational change and the accompanied electrogenic active transport are facilitated by the phosphorylation and dephosphorylation of an aspartate residue in the cytoplasmic domain, which is a hallmark of the P-type ATPase family (Axelsen et al., 1998).
The presence of the Na + /K + -pump is essential for excitability and secondary active transport. More than 40% of the energy produced in mammals is consumed by the Na + /K + -pump (Milligan et al., 1985). Although it is a machine designed for such a precise and important function, it has been shown that many cations, including congeners of K + and some organic cations, can bind to the same sites used by the pump to bind and transport K + ions (Mahmmoud et al., 2015;Ratheal et al., 2010). Competitive binding between Na + and other cations at the cytoplasmic side of the membrane has also been observed (Schneeberger et al., 2001). An unsolved puzzle, therefore, is how the Na + /K + -pump is able to recognize and accept two K + from the extracellular matrix, where Na + concentration is much higher, and how it selectively binds Na + from the cytoplasm to keep the pump cycle running forward.
Structural studies have provided tremendous insights into the function of the Na + /K + -pump, which consists of two obligatory subunits, a (catalytic) and b (auxiliary), and sometimes a tissue specific regulatory subunit from the FXYD protein family (Geering, 2006;Mercer et al., 1993). The transmembrane region of the a-subunit contains the ion binding sites within its 10 helices (called M1-M10). Recent crystal structures of the pump in its E 1 and E 2 states reveal the locations of the three Na + binding sites in E 1 and the two K + binding sites in E 2 (Kanai et al., 2013;Laursen et al., 2013;Morth et al., 2007;Shinoda et al., 2009). From structural alignments based on the heavy atom positions in the binding sites, it becomes clear that the binding pocket harboring sites I and II in E 1 overlaps with those in E 2 (Figure 1). Site III is only formed in E 1 . It is located between the transmembrane helices M5, M6, and M8 and is thought to exclusively bind Na + but appears to catalyze H + transport (Poulsen et al., 2010;Ratheal et al., 2010) in a manner that presents a complex dependence on the concentrations of Na + , K + , and H + (Mitchell et al., 2014). Previous studies have indicated that protonation at the ion binding pocket may play a role in the selectivity of these sites for K + when the pump is in its E 2 state (Yu et al., 2011). Biochemical assays on the mutant D926N, which is often used as a surrogate for protonated D926, also show that it induces distinct effects on Na + and K + binding (Einholm et al., 2010;Jewell-Motz et al., 1993), suggesting a potential change in its protonation state upon the E 1 /E 2 transition. Taken together, the evidence, although indirect, suggests the possibility of an E 1 specific protonation state that favors Na + binding to the pump. The nature of this protonation state is, however, unclear. eLife digest A protein called the sodium-potassium pump resides in the membrane that surrounds living cells. The role of this protein is to 'pump' sodium and potassium ions across the membrane to help restore their concentration inside and outside of the cell. About 25% of the body's energy is used to keep the pump going, rising to nearly 70% in the brain. Problems that affect the pump have been linked to several disorders, including heart, kidney and metabolic diseases, as well as severe neurological conditions.
The sodium-potassium pump must be able to effectively pick out the correct ions to transport from a mixture of many different types of ions. However, it was not clear how the pump succeeds in doing this efficiently. Rui et al. have now used a computational method called molecular dynamics simulations to model how the sodium-potassium pump transports the desired ions across the cell membrane.
The pump works via a so-called 'alternating-access' mechanism, repeatedly transitioning between inward-facing and outward-facing conformations. In each cycle, it binds three sodium ions from the cell's interior and exports them to the outside. Then, the pump binds to two potassium ions from outside the cell and imports them inside. Although the bound sodium and potassium ions interact with similar binding sites in the pump, the pump sometimes preferentially binds sodium, and sometimes potassium. The study performed by Rui et al. shows that this preference is driven by how protons (hydrogen ions) bind to the amino acids that make up the binding site.
The simulations also suggest that the pump uses a 'self-correcting' mechanism to prevent the pump from transporting the wrong types of ions. When incorrect ions are present at the binding sites, the pump cycle pauses temporarily until the ions detach from the pump. Only when the correct ions are bound will the pump cycle continue again.
In the future, Rui et al. hope to use long time-scale molecular dynamics simulations to show the conformational transition in action. In addition, the 'self-correcting' mechanism will be directly tested by letting the wrong and correct ions compete for the binding sites to see whether the pump will transport only the correct ions.
In the present study, we started with the recently published crystal structure of the Na + /K + -pump in a partially occluded Na 3 ÁE 1 (ADPÁPi) state and used molecular dynamics (MD) simulations to show that there is a correlation between the binding pocket protonation and the Na + selectivity in E 1 . The binding sites in Na 3 ÁE 1 (ADPÁPi) were tested one by one by 'alchemically transforming' the bound Na + into a K + while keeping the other two sites occupied by Na + . A similar approach was previously used to study the selectivity in other conformational states of the pump (Yu et al., 2011) (see also Materials and methods section). The results show that the Na + selectivity at all three sites is realized only when a specific set of binding pocket residues are protonated. This set of residues is, by comparison, different from that in the K + selective E 2 state. The implication is that the protonation state and the selectivity of the pump are tightly coupled; when the pump undergoes a transition between E 1 and E 2 , a protonation state switch occurs. The present findings also show that the effective selectivity of the pump is reinforced by a self-correcting mechanism, which prevents the occlusion step on either side of the membrane if the wrong type of ions were loaded into the binding sites. This mechanism ensures that counterproductive transport cycles do not occur.

pKa of the binding pocket protonatable residues
The crystal structures of the pump in its E 1 (PDBID 3WGV [Kanai et al., 2013]) and E 2 (PDBID 2ZXE [Shinoda et al., 2009]) conformations show that there is a large structural overlap between the Na + and K + binding pockets. The sites I and II are in the main binding chamber formed by helices M4, M5, and M6, while site III, located in between M5, M6, and M8 (Figure 1), is Na + exclusive and appears in E 1 only. The coordination of Na + and K + at sites I and II are similar. Many residues are found to coordinate both Na + and K + at these two sites in E 1 and E 2 ( Table 1). Although the structural difference between the ion binding pockets in E 1 and E 2 may not be particularly large, the local physicochemical environment can display considrable variations. In particular, the latter could affect the pKa and the protonation states of the six protonatable residues in the binding pocket ( Figure 1). Their pKa values calculated using PROPKA3.1 (Olsson et al., 2011) are listed in Table 2.
It is important to note that the pKa values calculated with empirical methods like PROPKA are sensitive to local structural perturbations. Even when the structures assume the identical conformational state and are taken from the same asymmetric unit from a crystal, the pKa values of the same residue can differ by more than one pH unit. For example, the pKa of E779 is 9.9 in chain A of the crystal structure 3WGV but the value is 8.4 in chain C, yet, the structural difference between the two chains is minimal ( Table 2). Because of this, only the structures with the highest resolution for the E 1 and E 2 states of the pump, 3WGV and 2ZXE, were used to guide the protonation state assignments in order to avoid false assignments of protonation state originated from structural uncertainty in lower resolution crystal structures. Interestingly, the protonation states of E327 and D804 assigned based on the 3WGV and 4HQJ structures are reversed. According to the PROPKA results, E327 appears to be deprotonated and D804 protonated in 3WGV, and vice versa in 4HQJ. While this inconsistency could be due to the lower resolution in the crystal structure of 4HQJ, it is worth noting that these two residues are in close proximity from one another suggesting that a proton could be passed back and forth between them in the E 1 state.

Protonation states and E 1 binding pocket stability
Previous calculations have indicated that a specific set of residues has to be protonated for the E 2 binding sites to be K + selective (Yu et al., 2011). Although the protonation states of D926 and E954 were not considered in the previous study because they lie outside of the two K + binding sites, they are likely to be protonated in E 2 according to the predicted pKa ( Table 2). The binding pocket  protonation derived from the crystal structure of the pump in its partially occluded Na 3 ÁE 1 (ADPÁPi) state differs from that in E 2 . According to their pKa values, D804, D808, and D926 should be deprotonated and the glutamates, E327, E786, and E954, should be protonated. Arguments based on pKa values predicted by PROPKA, however, have to be taken with caution, because the empirical method is highly sensitive to the local structural variation. For example, a deprotonated acidic side chain could have a PROPKA predicted pKa value at slightly higher than 7 because of the structural snapshot used to make the prediction. To seek a more robust assessment of these factors, molecular dynamics (MD) simulations were conducted to examine the structural stability of the binding pocket for different protonation configurations. One goal is to determine the protonation states of D804 and D926, both of which have a predicted pKa value within 1.5 pH unit to 7. The protonation state of D808 was also investigated since it affects the K + selectivity of the binding sites in the E 2 (K 2 ) state. A total of eight MD simulations were carried out (Table 3), starting from all possible protonation state configurations accessible by these three residues. The protocol for setting up these systems is given in details in the Materials and methods section. Figure 2 shows the binding pocket conformation in the E 1 systems ( Table 3) at the end of the allatom MD simulations. Out of the three aspartates when two or more are protonated the binding pocket becomes unstable, showing a large deviation from the crystal structure either with one of the bound Na + expelled from its binding site (systems E1_S5-7), or with a Clion attracted into the binding pocket (system E1_S4). These scenarios are not likely to happen in the normal function cycle of the pump. On the other hand, when the number of protonated aspartates is less than two the binding pocket remains structurally close to that in the crystal structure (systems E1_S0-3). The heavy atom root mean squared deviations (RMSD) between the crystal structure and the structure snapshots from the last 50 ns of the simulations are less than 1.6 Å . Table 3. Summary of the all-atom simulation systems and the FEP/H-REMD reduced systems. The binding site residues E327, E779, and E954 were kept protonated in all the systems.

Systems
Binding site residues Simulation time (ns) Protonation state and the E 1 binding pocket selectivity Analysis of the structural stability of the binding sites based on the MD simulations indicates that all of the four protonation states producing stable ion binding pockets may coexist when the pump is in the Na 3 ÁE 1 (ADPÁPi) state. The relative population of these states in the state Na 3 ÁE 1 (ADPÁPi), however, would differ from one another. The overall selectivity of the binding sites has contributions from all the states, and the predominant protonation state configuration should produce binding sites that are Na + selective. The determination of this protonation state configuration starts from the protein-membrane systems generated by the restraint-free MD simulations. A reduced structural model of the binding pocket is derived from each of these systems (Table 3) and the binding free energy differences between Na + and K + (DDG Na!K ) at the binding sites are calculated (See Materials and methods). The value of DDG Na!K reflects the affinity difference between K + and Na + binding as

MD
The results are plotted as the logarithm of the affinity ratio, ln K D;K =K D;Na À Á , in Figure 3. The values of DDG Na!K are shown in Table 4. As shown in Figure 3, the protonation state of the aspartates greatly affects both the E 1 and E 2 binding sites selectivity. The protonation state at E 2 (K 2 ) binding pocket proposed in the earlier study (Yu et al., 2011) is confirmed by the present calculations, in which the reduced system includes a larger region of the protein and more residues at the outskirt of the main binding pocket, including D926 and E954. It is evident that a different set of residues is protonated in the Na + selective E 1 as compared to in the K + selective E 2 ( Figure 3). While the glutamates (i.e., E327, E779, and E954) remain protonated in both the Na 3 ÁE 1 (ADPÁPi) and E 2 (K 2 ) states, the binding pocket aspartates take on opposite protonation states. In Na 3 ÁE 1 (ADPÁPi) D804 is protonated and both D808 and D926 are deprotonated, and their protonation states are reversed in E 2 (K 2 ). Free energy perturbation (FEP) calculations on the protonation/deprotonation of D804 estimated a pKa value of 9.2 (See Materials and methods), further supporting its protonated form under physiological conditions in E 1 .

Charge-neutralizing mutations and shift in selectivity
The impact of charge-neutralizing mutations of the key aspartate residues on the binding site selectivity was examined. D to N mutations have previously been used in experimental studies as a strategy to ascertain the possible charged state of these residues. The results of the free energy calculations, summarized in Figure 4, indicate that all the three sites remain Na + selective in the Figure 2. Comparison of snapshots at the end of the MD simulations (green) and the crystal structure of Na 3 ÁE 1 Á(ADPÁPi) (PDBID 3WGV) (white). The binding site residues are shown in stick presentation and the ions are shown as spheres. Binding site Na + ions from the MD simulation snapshots are in yellow, and the crystal Na + are in orange. A Clion has entered the binding site in system E1_S4 during the simulation and is shown in green. The view is from the extracellular side towards the intracellular side. The figure is produced with PyMOL (DeLano, 2002). DOI: 10.7554/eLife.16616.007 partly occluded Na 3 ÁE 1 (ADPÁPi) in both D804N and D808N, whilst D926N causes the sites to lose most of the Na + selectivity. Using Na + dependent phosphorylation and ATP binding assays, it was shown that the cytoplasmic binding affinity for both Na + and K + is reduced in D804N and D808N (Jorgensen et al., 2001;Pedersen et al., 1997). The D804N mutation affects the cytoplasmic K + and Na + binding differently. The mutation's impact on cytoplasmic Na + binding is not as prominent as on K + binding and D804N would appear more Na + selective than the wildtype. This is reflected in the calculation results. The K D,K /K D,Na at sites I and III are reduced by a factor of~50 from 476.6 and 2523.3 in the wildtype to 10.3 and 54.6 in the D804N mutant, but the sites are still Na + selective. The selectivity of site II is increased to a much larger extent, from 17.0 in the wildtype to 13359.7 in the mutant. This is an~1000 fold increase in selectivity and it implies that D804N is likely more Na + selective than the wildtype pump. The D808N mutant also becomes more Na + selective compared to the wildtype pump. The most dramatic increase in Na + selectivity happens at site II as in the case of D804N. Interestingly, site I, which prefers to bind K + in the wildtype with the protonated D808 ( Figure 3), is now Na + selective in D808N. It could be that the spatial packing is the major contributing factor to the site I ion selectivity. In the calculations, both sites II and III lose their Na + selectivity in the D926N mutant. The shifting of selectivity towards K + at site III in this mutant is a bit surprising, and worth commenting on. There are two possible explanations. First, substituting a negatively charged carbonyl oxygen with a bulkier but neutral -NH 2 at the D926 sidechain could prevent entrance of an ion to this site. Thus, the simulated conformation with an ion at this site could be energetically inaccessible. In other words, this is a metastable state with the absolute free energy of Na + and K + binding to this conformation equally prohibitive. An alternative explanation is that the D926N mutation alters the available conformational space accessible by helix M5 and this allows K + to go into site III as suggested in reference (Kanai et al., 2013). The binding of this K + , however, Figure 3. Ion binding sites selectivity characterized by ln K D;K =K D;Na À Á in states Na 3 E 1 Á(ADPÁPi) (blue) and E 2 (K 2 ) (red). Sites I (square), II (circle), and III (triangle) are distinguished by their shapes. Values from the previous calculations with a smaller reduced region are shown as empty symbols. All the binding site glutamates (i.e., E327, E779, and E954) are kept protonated. The protonation states of the binding site aspartates are indicated below the plot. DOI: 10.7554/eLife.16616.008 prevents the further binding of Na + , and results in the compromised Na + selectivity in this mutant. Even though the calculations show D926N with decreased Na + selectivity, they do not explain why experimentally the D926N mutant has compromised Na + binding ability in the absence of K + (Einholm et al., 2010). The loss in selectivity seen in this D926N mutant, however, could account for the strong inhibition by high K + on Na + /K + -ATPase activity (cf. Figure 2A in reference [Einholm et al., 2010]).
All the D to N mutations tested for the occluded E 2 (K 2 ) state appear to have little impact on the K + selectivity, in apparent contradiction with the experimental observations showing that both D804N and D808N display decreased selectivity for external K + (Kuntzweiler et al., 1996;Pedersen et al., 1997). In a previous computational study, the D808N mutation was also shown to minimally affect the K + selectivity in E 2 (K 2 ). Discrepancies between the calculations and the selectivity inferred from biochemical experiments have been noted previously, though the reason was not clear. A plausible explanation could be that the experimentally measured selectivity is an outcome from the relevant states including both P-E 2 and E 2 (K 2 ). However, any contribution from the P-E 2 state to the observed selectivity was not taken into account by the calculations because of the missing crystal structure of the outward-facing state. To test this idea, we generated a model structure of the P-E 2 state based on the homolog structures from the Ca 2+ SERCA pump (see below). This model was then used to calculate the DDG Na!K in the wildtype and the mutant pumps. The results show that both sites I and II have increased preference for external Na + in the P-E 2 state of D808N (Figure 4), which explains previous discrepancies between the calculations and the experiments.

Ion selectivity along the pump cycle
Using the crystal structures of the Ca 2+ SERCA pump as templates and a coarse grained transition pathway calculation method, ANMpathway (Das et al., 2014), we generated structural models of the Na + -and K + -loaded outward-facing P-E 2 state (see Materials and methods). The models appear to be structurally similar to the recently published P-E 2 like structures of the Na + /K + -pump (Gregersen et al., 2016;Laursen et al., 2013) and the vanadate-Inhibited, P-E 2 mimic of the Ca 2+ Table 4. The binding free energy difference (DDG Na!K ) at all the binding sites calculated from FEP/H-REMD simulations. The energy values are in kcal/mol.

Systems
Binding sites SERCA pump (Clausen et al., 2016). The MD simulations of the models are also able to predict with considerable accuracy the experimentally measured gating charge upon ion binding (Castillo et al., 2015). The structural transition pathways leading to these models provide a view of the intermediate structures along the pump cycle. Using these models, we calculated the DDG Na!K for the binding site of the intermediate state structures. The entire atomistic protein-membrane systems were used in the calculations. The results are presented in Figure 5 and Table 5. The calculations are valid as the DDG Na!K values mirror those computed from the reduced systems with the same protonation states ( Figure 3).
The results reveal a fascinating feature of the binding site selectivity along the pump cycle. Initially, the two K + binding sites in the outward-facing model P-E 2 ÁK 2 (state in Figure 5) appear to be non-selective or only weakly K + selective. However, the selectivity for K + over Na + increases as the pump occludes to form the K + -bound occluded state E 2 (K 2 ) (state ). After the dephosphorylation of P-E 2 , the binding pocket opens up to the cytoplasmic side. Accompanying this structural transition is the selectivity reversal at site I, switching the site from K + to Na + selective ( to in Figure 5). The protonation state change further shifts the ion selectivity at the sites in E 1 towards Na + . Changing the protonation state in K 2 E 1 ( to in Figure 5) to the protonation state dominant in Na 3 ÁE 1 (ADPÁPi), i.e., state in Figure 5 with protonated D804 and deprotonated D808 and D926, further reduces the K + selectivity at both sites I and II. When the pump is in this state, site I is Na + selective and site II is only weakly K + selective. Together, the E 2 to E 1 structural transition and the protonation state change promote the release of K + to the cytoplasmic side. Transient changes in the binding pocket protonation state upon occlusion/deocclusion are possible and could lead to variations in the magnitude of DDG Na!K for the intermediate states. Nonetheless, the general trend in free energy changes upon occlusion/deocclusion should remain because the DDG Na!K of the occluded state E 2 (K 2 ) has a much larger magnitude than that of the open-access outward-and inward-facing states.
The free energy calculations corroborate the notion that the two K + binding from the extracellular side is sequential and possibly cooperative. The sequential binding of extracellular K + was initially demonstrated by Forbush (Forbush, 1987). Recently using crystallography and isotopic measurements Ogawa et al. presented strong evidence that the first K + binds at site I and the second K + at site II (Ogawa et al., 2015). Whether there is cooperativity upon extracellular K + binding, however, is not clear from the crystal structures. Two electrode voltage clamp experiments have shown that the two extracellular K + binding events are relatively independent in the absence of extracellular Na + , but there is positive cooperativity of K + binding when extracellular Na + ions are present (Jaisser et al., 1994). The current DDG Na!K calculations of the two ion bound P-E 2 ( Table 5) provide an interpretation of this phenomenon. With a K + ion at site I, site II is relatively nondiscriminatory (DDG Na!K = 0.7 kcal/mol, p/p to p/s in state P-E 2 .K 2 , Table 5), but it becomes much more Na + selective (DDG Na!K = 2.4 kcal/mol, s/s to s/p in state P-E 2 ÁNa 3 , Table 5) when there is a Na + ion bound at site I. In the presence of extracellular Na + , K + ions have to first compete with Na + to bind at site I in order for the subsequent K + to bind, therefore resulting in the observed binding cooperativity in the presence of extracellular Na + . A similar phenomenon is seen in the Na + branch of the pump cycle. Initially, site I in the inwardfacing Na 3 ÁE 1 (ADPÁPi) structure (state ) is weakly Na + selective when the other two sites are filled with Na + . However, the intermediate state during the occlusion process shows that its site I is highly discriminating against K + with a DDG Na!K of 2.9 kcal/mol ( Figure 5 and Table 5). Therefore, if a K + is bound in this site, the free energy barrier toward occlusion would be much higher than when a Na + ion was bound and the subsequent occlusion is not likely. In the case when a K + ion replaces a Na + at site II or III, although the energy barrier of occlusion is not as prominent, such a binding pocket ion configuration is energetically less favorable than the native configuration and its appearance is unlikely in the first place.
Unlike the K + branch, the DDG Na!K calculations along the Na + branch offers limited insights on the sequence of Na + binding from the cytoplasmic side. The results are indicative of the selectivity at the sites, but not the absolute binding affinity. Based on the calculations alone, it is not possible to determine which site is the first to bind cytoplasmic Na + . It could be site III, as suggested by Kanai et al., and this prepares the other two sites for the following Na + binding (Kanai et al., 2013). Or, alternatively, the first two Na + bind at sites I and II in the main binding pocket and the last Na + enters, takes over site II, and pushes the two previously bound ions further, so that the ions in sites II and I now move to sites I and III, in a process reminiscent of the "knock-on" mechanism occurring in potassium channels (Hodgkin et al., 1955). A direct assessment of these proposed binding sequences will require further experiments.
The heavy atom RMSD between the binding pocket residues is 2.5 Å and the structural difference remains after hundreds of ns of simulations. Empirical pKa and FEP calculations based on the MD simulation equilibrated structures indicate that the binding pocket glutamates (i.e., E327, E779, and E954) are likely protonated in both Na 3 .E 1 Á(ADPÁPi) and E 2 (K 2 ), although previous mutagenesis experiments showed that charge-neutralizing mutation E327Q affects pump function, possibly by altering ion-pump interactions and the kinetics of the occlusion/deocclusion reactions along the pump cycle (Jorgensen et al., 2001;Kuntzweiler et al., 1995;Li et al., 2006;Nielsen et al., 1998). The calculations also suggest that the protonation states of D804, D808, and D926 are different in Na 3 ÁE 1 Á(ADPÁPi) and E 2 (K 2 ). Among the three aspartates only D804 is protonated in order for all three sites to stay Na + selective in Na 3 ÁE 1 Á(ADPÁPi). The E 1 binding pocket devoid of ions carries a net charge of -2, consistent with that deduced from previous fluorescence studies (Schneeberger et al., 1999). Mutagenesis experiments are also consistent with the unlikely protonation of D808 and D926 when the pump is in E 1 trying to bind cytoplasmic Na + (Jewell-Motz et al., 1993;Pedersen et al., 1997). The protonation states of these three residues are reversed in E 2 (K 2 ) with D804 deprotonated and the other two protonated, a result that is supported by a previous computational study (Yu et al., 2011). The different protonation states for E 1 and E 2 also agree well with the 'four-site' model proposed by Skou and Esmann more than 30 years ago with the K + -bound E 2 state carrying two sidechain protons (H 2 EK 2 ) and the Na + -bound E 1 state carrying only one (HENa 3 ) (Skou et al., 1980). It is worth pointing out that the second proton (i.e., the proton on D926) in the K + bound E 2 state must come from and return to the same side of the membrane during the pump cycle so that the net charge moved per cycle by the pump is one. It seems unlikely that this proton could come from the extracellular side because altering extracellular pH would then protonate or deprotonate D926, causing major changes in Na + binding affinity and the maximum pumping turnover rate. This is not observed over an external pH range 9.6 to 5.6 (Mitchell et al., 2014;Vasilyev et al., 2004;Yu et al., 2011). Therefore, it is more likely that the D926 proton comes from the cytoplasmic side. This is supported by MD simulations of the P-E 2 ÁNa 2 and Na 2 ÁE 1 Á(ADPÁPi) revealing the existence of aqueous pathways connecting the cytoplasm and D926 (Figure 6). One of the water pathways is located between the helices M5, M7, and M8 ( Figure 6B), similar to the C-terminal proton pathway previously proposed (Poulsen et al., 2010). A proton could enter through this pathway to protonate D926 and then leave through the same pathway during the E 2 to E 1 transition, or through an alternative path passing the main binding pocket along with the dissociating ions as seen in the Na 2 ÁE 1 Á(ADPÁPi) simulation ( Figure 6A).
Under physiological conditions, the challenge faced by the Na + /K + -pump is to effectively pick out the correct ion species from a solution much more concentrated with other types of ions. Even when the binding affinity of the ion species being transported is a few orders of magnitude higher than that of the other ion (DDG bind = 1 to 3 kcal/mol), the advantage in binding is undermined by the higher concentration of the competing ion. The pump overcomes this problem by raising the free energy barrier for the occlusion step in the presence of incorrect ions in the binding pocket, thus preventing the faulty transport ( Figure 5 and Table 5). This self-correcting mechanism makes sure that the intracellular and the extracellular gates do not close (i.e., occlude), which is required for the pump cycle to go forward, unless the correct ion configuration is present at the binding pocket. Therefore, this is an inherent part of the ping-pong mechanism the pump uses to transport ions with high selectivity. Although other energetic barriers, like pathways the ions use to travel to the binding site, could contribute to the ion binding specificity, it is not likely the case here as both Na + and K + can access the binding pocket in the simulations of the E 1 and P-E 2 states without ions bound.
The concept of a self-correcting pumping mechanism sheds new light on a number of puzzling functional observations. For instance, it explains why E 1 displays a strong apparent affinity for K + in competitive binding assays (Schneeberger et al., 2001), even though the pump still binds and occludes Na + from the cytoplasmic side. The calculations show that all the sites in the Na 3 ÁE 1 Á(ADPÁPi) state are selective for Na + , likely because this is a partially occluded state. A wrong combination of binding site ions would not have been occluded and reached such a state. In a fully inward-open conformation, the selectivity of the sites ought to be fairly weak. The high affinity for K + observed in E 1 in these experiments is due to the backward occlusion leading to the K + -bound state E 2 (K 2 ). The proposed self-correcting mechanism parallels the suggestion based on experiments that the phosphorylation and occlusion of E 1 requires 3 Na + bound, and this increases the apparent affinity for Na + in the normal pump cycle (Schneeberger et al., 2001). The self-correcting mechanism can also account for the different effects of Na + and K + congeners on the release rate of the occluded 86 Rb + (a K + congener) to the extracellular side upon the "backdoor" phosphorylation in the presence of Pi reported by Forbush (Forbush, 1987). This backward deocclusion is accelerated in the presence of Na + , and it follows a single phase, while a much slower deocclusion process with a second slow phase is observed in the presence of K + or Rb + . According to the self-correcting mechanism, it is likely that upon the replacement of one of the two occluded 42 K + or 86 Rb + by cold K + at one site, the pump succeeded to occlude again, resulting in the second slow phase. These explanations are sensible, in turn lending support to the proposed self-correcting mechanism. Given that the selectivity in the outward/inward facing states is after all not as strong as previously thought, such a mechanism must be in place. This analysis suggests that an important component of the overall selectivity emerges from the increased free energy barrier associated with the occlusion process while the wrong type of ion is bound to a site. Although the pump would be kinetically more efficient by preventing the wrong ions from reaching the binding sites to begin with, the free energy calculations do not support the notion of highly selective binding sites for the open-access states. While such mechanism may seem inefficient because the pump must try to figure out that there is an issue with selectivity only after binding of several ions of the wrong type, it is important to realize that the Na + /K + -pump is not a particular fast molecular machine. The estimated turnover rate is less than 100 per second and decreases even further as the transmembrane voltage becomes more negative (Heyse et al., 1994). The implication is that the system has plenty of time, at the molecular level, to function near thermodynamic equilibrium. In other words, the pump has not been evolutionarily optimized to be a particularly fast molecular machine, but an energetically efficient one. Thus, even though the idea of a self-correcting ion selectivity mechanism seems counterintuitive and inefficient, it is consistent with the physical conditions under which the pump has to operate.
In summary, the present study highlights the tight coupling between the Na + /K + selectivity of the binding sites, the protonation state of the coordinating residues, and the conformational state of the pump. The important functional consequence of such tight coupling is the necessity to have the correct type and number of ions in the binding pocket for the pumping cycle to proceed forward toward the occlusion step. Because the ion binding selectivity is strongly dependent upon the protein conformation, a self-correcting mechanism counteracting the effect of the ion concentration in the environment ensures an efficient function of the pump.

Materials and methods
Preparing the all-atom membrane-Na + /K + -pump simulation systems The crystal structures, 3WGV (Kanai et al., 2013) and 2ZXE (Shinoda et al., 2009), representing the Na + /K + -pump in its Na 3 ÁE 1 (ADPÁPi) and E 2 (K 2 ) states were used to build the simulation systems. The structure 3WGV contains two copies of the pump assembly in the asymmetric unit. Since the structural variations between the two assemblies are minimal, only the copy including chains A, B, and G was kept. Several small molecule ligands were co-crystalized with the pump in 3WGV, including an ADP, a AlF À 4 ion, an oligomycin A, two Mg 2+ , three cholesterol, four Na + , and five 1,2-diacyl-sn-glycero-3-phosphocholine molecules for each copy of the pump structure. Among these ligands, oligomycin A was not included in the simulation system. The ion was replaced by a PO 3À 4 and POPC molecules were used in place of the 1,2-diacyl-sn-glycero-3-phosphocholine. The structure of 2ZXE also contains small molecule ligands including cholesterol, K + , Mg 2+ , and MgF 2À 4 . Similarly, a PO 3À 4 was used to replace the MgF 2À 4 ion. Since the Na + /K + -pump crystal structures were obtained from different organisms, the residue numbers differ slightly. For the sake of convenience, the numbering scheme in the newly resolved 3WGV from pig kidney was adopted in this study.
After removing irregularities from the PDB files, the Na + /K + -pump subunits were capped with acetylated N-termini and amidated C-termini. The ectodomain of the b-subunit was not included in simulation systems to reduce the computational cost. The orientation was chosen to be the same as in the OPM database (Lomize et al., 2006). At this stage, the protonation states of the binding site residues were assigned with the PATCH command in CHARMM (Brooks et al., 2009). Eight different protonation states were considered in Na 3 ÁE 1 (ADPÁPi) state systems and both the protonated and charged D926 were included in the E 2 (K 2 ) state systems. This resulted in a total of ten systems before proceeding to the next step. Table 3 shows the system names and the associated binding site protonation. We used the membrane builder module in CHARMM-GUI (Jo et al., 2007(Jo et al., , 2008(Jo et al., , 2009Wu et al., 2014) to generate POPC bilayers around the pump structures. Once this was completed, each protein-membrane complex was solvated by an equal molar mixture of KCl and NaCl. The final system had a combined cation concentration of 0.3 M (i.e., [K + ] = [Na + ] = 0.15 M). At the end the E1 system was 84 Â 110 Â 158 Å 3 in size and contained~138,000 atoms, while the dimension of the E2 system was 85 Â 109 Â 105 Å 3 with~134,000 atoms. Each completed system was subjected to a 675-ps equilibration with reducing restraints on the heavy atoms to relax the initially uncorrelated components, followed by a 10-ns unrestrained pre-production using the simulation package NAMD2.9 (Phillips et al., 2005). After the systems were well equilibrated, they were simulated longer using the special-purpose supercomputer Anton (Shaw et al., 2009), which is designed for long time scale MD simulations.

Generating the D to N mutant systems
Experimentally asparagine and glutamine are used as surrogates for protonated aspartate to study the effect of protonation. The effects of these charge-neutralizing mutants on the selectivity of the pump can be investigated computationally. The most interesting mutations in the context of this study are the single D to N mutations at the binding pocket, including residues D804, D808, and D926. The mutations were made on the crystal structures by replacing the proton on the OD2 atom in the aspartate with an -NH 2 amine group. D804N and D808N mutations were also made in the outward facing P-E 2 structural model (see below) to study how they affect the external K + binding. The protonation states of the other titratable residues were determined with PROPKA. Each of these systems was equilibrated without any restraints for 40 ns. The system snapshot with the least structural deviation of the pump to the averaged structure during the simulation was used to generate the reduced system. The mutant systems are shown in Table 3.

Selectivity at the binding site
The absolute free energy of an ion i binding to a binding site inside a protein has the following form, The difference in the first term, G site i; int À G bulk i; int , represents the nonbonded interaction component of the free energy change upon moving the ion from the bulk solution to the binding site. The subtraction in the second term, Àk B TIn F t C ð Þ À G site i; trans , reflects the lost of translational freedom. The translational freedom factor F t in bulk solution can be evaluated analytically under a rigid rotor approximation (Deng et al., 2006). Its final form depends on the force constants and the equilibrium values in the distance and angle restraints applied on the ion and the surrounding protein atoms. Based on Equation (1), the selectivity of a binding site can be expressed as the binding free energy difference of two ion species. For example, the binding free energy change upon changing ion i to j at the binding site is A negative value of DDG i!j indicates that ion j binds more favorably than ion i and the site is j selective. There are three terms to be evaluated in Equation (2). DG site i!j and DG trans i!j are calculated using the reduce binding site model, while DG bulk i!j is computed using a water sphere with the impact from the bulk solution factored in with a boundary potential (see below).

Generating the reduced binding site systems
To generate the reduced binding site, the all-atom system prepared according to the procedures above was divided into an inner region and an outer region. The inner region was defined as residues and water molecules within 15 Å to the center of mass of the bound ions. Everything within this region was treated explicitly. An extended inner region was specified by a 3-Å thick shell continuing from the inner region outwards to create a smooth spherical dielectric cavity. Water molecules in this region were removed and their impact on the binding site was included implicitly. The coordinates of atoms in these extended region and those linked to them within three atomic bonds in the inner region were held fixed during the simulations. The rest of the system was considered the outer region, in which the atoms were removed and their impact on the inner region atoms was represented by the General Solvent Boundary Potential (GSBP) in the form of a solvent-shielded static field and a solvent-induced reaction field (Im et al., 2001). The reaction field arising from changes in charge distribution in the inner region was expressed in terms of a generalized multipole expansion using 11 spherical harmonic functions. Both the solvent-shielded static field and the reaction field matrix were independent of the inner region configuration, and therefore were calculated only once before further simulations using the finite-difference PoissonÀBoltzmann (PB) method with the PBEQ module (Im et al., 1998) in CHARMM. In these calculations a dielectric constant of 1 was used for the inside of the protein within the inner and outer regions, whereas the rest of the system had a dielectric constant of 80. The atomic Born radii used in the PB calculations were determined by free energy calculations in explicit solvent (Nina et al., 1997). All these reduced systems were further equilibrated for 200 ps using Langevin dynamics at 303.15 K with a friction coefficient of 5 ps À1 assigned to all non-hydrogen atoms. All bonds involving hydrogen atoms were fixed with the SHAKE algorithm (Ryckaert et al., 1977). Nonbonded interactions within 14.5 Å were accounted for explicitly, while everything else beyond this distance was treated with an extended electrostatics method (Stote et al., 1991). The simulation program CHARMM (Brooks et al., 2009) was used for the equilibration. Table 3 lists all the reduced systems.

Free energy perturbation simulations with boosting potential
It is known that the binding free energy calculated from FEP/MD simulations suffers from convergence issues because the residue sidechains at the binding site only sample a few of the several accessible rotameric states. This problem can be alleviated by introducing a replica-exchange scheme aiming at enhancing the sampling of sidechains (Jiang et al., 2009(Jiang et al., , 2010. This scheme allows exchange between the neighboring l windows. Each l window is coupled with a set of replica systems with a boosting potential of increasing strength used to accelerate the inter-conversion between sidechain rotameric states. We employed this hybrid FEP/H-REMD scheme implemented in the REPDSTR module in CHARMM (Jiang et al., 2010) to calculate the DDG i!j .
The boosting potential for the c 1 sidechain torsion of a binding site residue was obtained by fitting the potential of mean force as a function of the torsion c, WðÞ, with a cosine Fourier series in the form of The angle c 1 is dihedral formed by the bonded atoms N, CA, CB, and CG. The total number of the cosine terms, N, varies from 3 to 6, depending on which one produce a better fit to the WðÞ. A residue is considered in the binding site if any of its sidechain heavy atoms is within 4.5 Å to the bound Na + or K + . Table 6 lists all the binding site residues included in the boosting potential calculations.
We performed umbrella sampling simulations to obtain the WðÞ for each binding site residue. First, the entire transmembrane helix containing the residue of interest was taken from the crystal structure with its orientation kept the same as in the OPM database. This was to provide a proper secondary structure environment. Next, we replaced the sidechains of all the other residues on the helix to hydrogen atoms to remove their steric effects. An implicit planar membrane model, EEF1/IMM1 (Lazaridis, 2003;Schneeberger et al., 1999) was used in place of a membrane bilayer to solvate the helix. 72 umbrella windows along c were set up with a window spacing of of 5˚. A harmonic bias potential was applied in each window with a force constant of 100 kcal/mol/rad 2 . After the windows were generated, each of them was simulated for 50 ns at 303.15 K using Langevin dynamics with a 2 fs time-step. All bond lengths involving hydrogen were fixed with the SHAKE algorithm (Ryckaert et al., 1977). The cutoff distance for nonbonded interactions was set to 11 Å . The simulation program CHARMM (Brooks et al., 2009) was used to conduct the simulations. The resulting distributions of c were unbiased using the weighted histogram analysis method (WHAM) (Kumar et al., 1992). Figure 7 shows the WðÞ and the fitted curves. The fitting parameters and are given in Table 6. An appropriate boosting potential can be easily applied using the CONS DIHE command in CHARMM, with the sign in front of U BP ðÞ reversed. This will effectively cancel out the potential barrier between the rotamers of a given residue sidechain.
Before computing DG site i!j and DG trans i!j , restraints were set up to restrict the translational freedom of the bound ion of interest. First, three points within the protein region were picked out. They were used along with the position of the bound ion to set up the restraints. We employed the same protocol as in the Ligand Binder module (Jo et al., 2013) in CHARMM-GUI to select the three protein anchor points, p 1 , p 2 , and p 3 . The relative position of the bound ions to the protein can be defined by the distance r from l 1 to p 1 , the angle between l 1 , p 1 and p 2 , and the torsion along l 1 , p 1 , p 2 , and p 3 . The translational restraint potential took the form of The force constant regarding distance (i.e., k r ) was set to 10 kcal/mol/Å 2 , and the rest of the force constants were 200 kcal/mol/rad 2 . The equilibrium values in u trans were taken from the 200 ps equilibration of the reduced binding site.

DG site
Na!K was computed with a set of FEP/H-REMD simulations. The translational restraint, u trans , was applied to restrict the translational freedom of the ion. The Na + ion was changed alchemically to a K + with an alchemical coupling factor, l. 16 l windows were used. Each l window included 8 replicas with the strength of boosting potential scaled from 0 to 1. Exchange attempts were made every 0.2 ps and were only allowed between neighboring replicas with different boosting potential strengths in the same l window and between the neighboring l windows with zero boosting potential. A total of 128 replica systems were included in the calculations. Each replica system was simulated for 2 ns. To compute DG trans Na!K at a given binding site, two sets of simulations were set up, one with Na + at the binding site and the other with K + for the calculation ofDG site Na; trans and DG site K; trans . The translational and orientational restraints were decoupled gradually (Beglov et al., 1994) using a coupling factor l. The l window and the boosting potential setup were similar to those used in the DG site Na!K calculations. All the FEP/H-REMD calculations were performed using the REPDSTR module in CHARMM (Jiang et al., 2010). The output energies from the zero boosting potential systems were collected and processed using WHAM (Kumar et al., 1992). To compute the standard error of DDG, we divided the trajectories into 10 blocks. The standard deviations of the averaged DDG from all the blocks are computed and reported in Tables 4 and 5 and Figures 3-5. Cautions should be taken when comparing the calculated and experimentally measured DDG, because the latter usually Table 6. Binding site residues and the fitted parameters k i and c 0,i . contains contributions from multiple states of the pump while the calculated DDG is done using a single state.
The bulk ion-water system and the calculation of DG bulk Na!K The bulk system was generated by building a water sphere of 10 Å radius and centering it at the origin. The water sphere was made form pre-equilibrated water boxes with TIP3P water molecules. A Na + ion was placed at the center of the sphere. Any water molecules overlapping with the ion were deleted. The influence of the remaining bulk was approximated by the spherical solvent boundary potential (Beglov et al., 1994). The system was equilibrated for 200 ps at 303.15 K. Other simulation options were kept the same as described in the reduced binding site system. During the equilibration the position of the ion was restrained using a weak harmonic bias potential with a force constant of 0.5 kcal/mol/Å 2 . Once equilibrated, the Na + was gradually changed to a K + using the PERT module in CHARMM with 11 l windows. The simulation time for each l window was 1 ns. The resulting DG bulk Na!K is 18.34 kcal/mol after unbiasing the energy outputs with WHAM. The calculated DG Na!K for all the binding sites in the systems summarized in Table 3 are shown in Table 4.

D804 pKa calculations with explicit solvent
To further confirm the protonation state of D804, we evaluated its pKa shift with additional simulations in explicit solvent using the following formula: are the free energy change of aspartate deprotonation in the protein environment and in bulk water, respectively. The reduced system of E1_S1 (Table 3) was used to compute the deprotonation free energy of D804 at the ion binding site (DG deprot site ). To calculate the DG deprot bulk , an aspartic acid residue with an acetylated N-terminus and an amidated C-terminus was put into a pre-equilibrated water sphere of 15 Å radius. As in the calculations of DG bulk Na!K and DG bulk K!Na , the impact of the bulk solution beyond the current system was incorporated by the spherical solvent boundary potential (Beglov et al., 1994). Alchemical FEP calculations were carried out in both systems. The l windows were evenly spaced to gradually deprotonate the aspartate. The numbers of windows used in the DG deprot site and DG deprot bulk calculations were 24 and 10 respectively. Each window was equilibrated for 200 ps and contained 5 ns of sampling. The calculated DG deprot site and DG deprot bulk are À44.8 kcal/mol and À51.9 kcal/mol, respectively. Using Equation (5), the DpKa is 5.1, meaning the pKa of D804 is right shifted by 5.1 pK unit when it is in the Na + /K + pump ion binding site. The final calculated pKa of D804 is 4.1 + 5.1 = 9.2 with 4.1 being the pKa value of an aspartate in bulk water (Berg et al., 2002), reinforcing the conclusion that D804 is protonated in the E 1 state of the pump.

Transition pathway calculations
Three conformational transition pathways were generated based on the shared homology between the Na + /K + -pump and the Ca 2+ SERCA pump. One of such pathways connects states E 2 (K 2 ) and P-E 2 ÁK 2 . The second connects states E 2 (K 2 ) and K 2 ÁE 1 , and the third describes the transition between states Na 3 ÁE 1 -P and P-E 2 ÁNa 3 . First, A Ca-atom-only transition pathway for each of these was generated using the ANMPathway online server (http://anmpathway.lcrc.anl.gov/anmpathway.cgi) (Das et al., 2014) based on the SERCA pump crystal structures including 3B9B (Olesen et al., 2007), 1WPG , and 1VFP . Since both the Na + /K +pump and the Ca 2+ SERCA pump are P-type ATPases, it is reasonable to assume that they have similar pumping mechanisms and they share the same set of states along the pump cycle. Although crystal structures of the Na + /K + -pump are scarce, multiple high-resolution structures of the SERCA pump in different states are available (Karlsen et al., 2016). Among them, 3B9B represents the outward facing P-E 2 state. All the other SERCA pump structures were aligned based on their transmembrane region Ca positions to the Na + /K + -pump and the Ca RMSD were computed to find those that resemble the structural states captured by the Na + /K + -pump crystal structures 3WGV (Kanai et al., 2013) and 2ZXE (Shinoda et al., 2009). The two SERCA pump structures representing states Na 3 ÁE 1 ATP (3WGV) and E 2 (K 2 ) (2ZXE) are 1VFP and 1WPG, respectively. Each generated coarsegrained pathway is made of a sequence of structural snapshots containing Ca atoms only. These snapshots are called images and are distributed at equal RMSD intervals. The images from the coarse-grained pathways were then used in the all-atom targeted molecular dynamics (TMD) simulations to guide each transition. The starting system in the first and second pathway TMD simulations was taken from the well-equilibrated MD simulation system E2_S1 starting from the crystal structure 2ZXE. The third pathway was realized using the MD simulation system E1_S1 built from the crystal structure 3WGV, with the catalytic D369 phosphorylated. The protocol of the TMD simulations followed those published before (Castillo et al., 2015) Selectivity calculations along the pump cycle Binding site ion selectivity was evaluated for nine systems representing different stages along the pump cycle. Among these systems are well-defined states: P-E 2 ÁK 2 , E 2 (K 2 ), K 2 ÁE 1 , Na 3 ÁE 1 (ADPÁPi), and P-E 2 ÁNa 3 . The structures used for them are the generated P-E 2 ÁK 2 model, the equilibrated crystal structure 2ZXE, the equilibrated crystal structure 3WGV with bound Na + and with bound Na + replaced by K + at sites I and II, and the generated P-E 2 ÁNa 3 model. Several intermediate states were also included. Intermediate state is between states P-E 2 ÁK 2 and E 2 (K 2 ), taken from the TMD simulations of the first path after image 35. State was taken at image 33 along the transition from state E 2 (K 2 ) to K 2 ÁE 1 . Along the same vein, state represented the intermediate at image 29 between states Na 3 ÁE 1 (ADP.Pi) and P-E 2 ÁNa 3 . The intermediates were taken near the midpoints of the transitions. An additional state K 2 ÁE 1 * with altered binding site residue protonation was also included to reflect the protonation state change upon the E 2 to E 1 transition. All the states included in the selectivity calculations and their relative placement along the pump cycle can be seen in Figure 5.
To evaluate the selectivity at a given binding site in a system, three simulations were performed with slightly varied Lennard-Jones (LJ) parameters for the ion of interest. In the first simulation, the parameters of the ion were unaltered. A linear interpolation of the LJ and the NBFIX parameters between Na + and K + was made. In the second simulation, the LJ and NBFIX parameters of the ion were replaced by those of a Na/K hybrid at the middle point of the interpolation. In the last simulation, the ion's parameters were completely changed to represent the other ion type, either K + or Na + , depending on the starting ion type. 10 ns trajectories were generated for each simulation using the simulation package OpenMM6.2 (Eastman et al., 2013). Energies were calculated from each trajectory using two parameter sets. The difference was fed into the WHAM equation (Kumar et al., 1992) and the free energy changes upon mutating the ion to the intermediate as well as upon mutating the intermediate to the other ion type were solved self-consistently. The sum of the two free energies gives the DG site . The same strategy was used to compute the DG bulk in water (DG bulk kcal/mol). The binding free energy difference is given by DDG bind ¼ DG site À G bulk . The selectivity of binding site with different occupancies was also included in the calculations. The DDG bind are shown in Table 5.

NAMD simulation protocol
The initial relaxation and the restraint-free equilibration of the membrane pump systems were performed using the NAMD2.9 (Phillips et al., 2005) simulation package with the input scripts from the CHARMM-GUI membrane builder module. The NAMD simulation temperature was set to 303.15 K using Langevin Dynamics with a damping coefficient of 10 ps À1 during the relaxation and 1 ps À1 during the restraint-free equilibration. The van der Waals interactions were smoothly switched off at 10-12 Å by a force switching function (Steinbach et al., 1994) and the electrostatic interactions were calculated using the particle-mesh Ewald method with a mesh size of~1 Å .

Anton simulation protocol
After a short equilibration with NAMD, each all-atom system listed in Table 3 were subjected to a hundred ns scale long production without any restraints using the special-purpose supercomputer Anton (Shaw et al., 2009). The volume of the periodic cell was kept constant and the temperature was set to 303.15K using the Nosé -Hoover thermostats . The lengths of all bonds involving hydrogen atoms were constrained using M-SHAKE (Kräutler et al., 2001). The cutoff of the van der Waals and short-range electrostatic interactions was set to an optimal value suggested by the Anton guesser script, guess_chem. Long-range electrostatic interactions were evaluated using the k-space Gaussian split Ewald method (Shan et al., 2005) with a 64 Â 64 Â 64 mesh. The integration time step was 2 fs. The r-RESPA integration method  was employed and long-range electrostatics were evaluated every 6 fs.

OpenMM simulation protocol
Simulations used to compute ion selectivity along the pump cycle were conducted using the simulation package OpenMM6.2 (Eastman et al., 2013). The constant pressure and temperature (NPT) ensemble was used for these simulations. The pressure was maintained at 1 atm with a Monte Carlo barostat, which attempts to adjust the system volume every 0.2 ps. The Langevin dynamics algorithm with a 1.0 ps À1 friction coefficient was used to hold the simulation temperature, which is at 303.15 K. The lengths of all bonds involving hydrogen were fixed and the integration time step was set to 2 fs. A force switching function was applied from 10 to 12 Å to gradually turn off the van der Waals interactions and the particle-mesh Ewald method with an error tolerance of 0.0005 was used to evaluate the electrostatic interactions.

Simulation force field parameters
The same force field parameters were used in all other simulations in this study. The PARAM27 allatom force field of CHARMM (MacKerell et al., 1998) with a modified version of dihedral crossterm correction (Mackerell et al., 2004) was used for the protein and the C36 lipid force field (Klauda et al., 2010) was used for POPC. Water molecules were modeled with the TIP3P potential (Jorgensen et al., 1983).