Molecular simulations and free-energy calculations suggest conformation-dependent anion binding to a cytoplasmic site as a mechanism for Na+/K+-ATPase ion selectivity

Na+/K+-ATPase transports Na+ and K+ ions across the cell membrane via an ion-binding site becoming alternatively accessible to the intra- and extracellular milieu by conformational transitions that confer marked changes in ion-binding stoichiometry and selectivity. To probe the mechanism of these changes, we used molecular simulation and free-energy perturbation approaches to identify probable protonation states of Na+- and K+-coordinating residues in E1P and E2P conformations of Na+/K+-ATPase. Analysis of these simulations revealed a molecular mechanism responsible for the change in protonation state: the conformation-dependent binding of an anion (a chloride ion in our simulations) to a previously unrecognized cytoplasmic site in the loop between transmembrane helices 8 and 9, which influences the electrostatic potential of the crucial Na+-coordinating residue Asp926. This mechanistic model is consistent with experimental observations and provides a molecular-level picture of how E1P to E2P enzyme conformational transitions are coupled to changes in ion-binding stoichiometry and selectivity.

Na ϩ /K ϩ -ATPase is a membrane protein that actively transports sodium ions out of the cell while importing potassium ions, both against their electrochemical gradients, thus providing potential energy necessary for many essential cellular functions (1)(2)(3). Malfunction of Na ϩ /K ϩ -ATPase has been linked to numerous diseases including impaired memory and learning, familial hemiplegic migraine 2, rapid-onset dystonia Parkinsonism, and heart failure (4 -6); thus, Na ϩ /K ϩ -ATPase is an important target for treatment of brain and heart conditions (7,8).
By harnessing chemical energy stored in ATP, the Na ϩ /K ϩ -ATPase cycles between two major conformational states during active ion pumping: one state with high affinity for sodium ions (E1) and a second with high affinity for potassium ion (E2). Although the complete mechanism of function of the Na ϩ /K ϩ -ATPase is more complicated, involving several conformational transitions described more fully in the Post-Albers scheme (9 -11), here we focus on the sodium-bound E1P and potassiumbound E2P states (i.e. phosphorylated E1 and E2) for which crystal structures are available (12)(13)(14)(15)(16)(17). One of the central features of ion transport by the Na ϩ /K ϩ -ATPase is a change in ion selectivity between E1 and E2 conformations (3). The origin of this selectivity change has been a focus of Na ϩ /K ϩ -ATPase research since its discovery by Jens Christian Skou in 1957 (2). Extensive mutational studies have identified key residues involved in ion binding, which include five acidic residues, Asp 804 , Asp 808 , and Asp 926 , Glu 327 , and Glu 779 (13), whose orientation, distance, and ion coordination have been proposed to be involved in determining selectivity (3). The first crystal structure of the Na ϩ /K ϩ -ATPase, solved for the potassiumbound E2P state in 2007 (15), and later followed by several higher resolution structures (14,16,17) revealed a binding site in which these five acidic residues and several backbone carbonyl oxygens tightly coordinate two potassium ions within an ϳ15-Å 3 space. The close proximity of acidic binding residues ( Fig. 1) suggested that some of these residues must be protonated. Indeed, protonation of binding site residues is intrinsic to the transport cycle of the Ca 2ϩ pump, sarco/endoplasmic reticulum Ca 2ϩ -ATPase, a closely related P-type ATPase (18).
Ion-binding selectivity and stoichiometry are central to the efficient countertransport of Na ϩ and K ϩ by the Na ϩ /K ϩ -ATPase, and such properties likely depend on the protonation state of acidic residues in the ion-binding site. The notion that protonation states can modulate ion selectivity is strongly supported by a joint computational and experimental study by Yu et al. (19) in which different protonation states of the potassium-bound E2P state were predicted to confer differences in ion selectivity. In this work, electrophysiological experiments con-firmed computational predictions that potassium selectivity decreases with increasing pH, consistent with previous pH studies (20). However, the exact protonation state could not be determined nor could the results provide insight into the protonation state of sodium-bound conformations. Most importantly, it is not clear how changes in protonation state are controlled by or coupled to structural changes during active transport of ions. The first crystal structures of the sodiumbound E1P state were obtained in 2013 (12,13), providing an unprecedented opportunity to discover the molecular origins of selectivity in atomic detail, now a half-century-old question.
Here, we address the question of how ion selectivity and binding stoichiometry are achieved at the ion-binding sites in different protein conformations by extensively scrutinizing Na ϩ /K ϩ -ATPase in atomic detail using molecular simulation. In particular, we use pairwise additive force fields (CHARMM36) within the framework of several molecular dynamics simulation-based approaches, e.g. metadynamics and free-energy perturbation. Our physical model of protein-ion interactions does not take into account the polarization of the ion's ligands; as such, the molecular description underlying our molecular simulations is only moderately accurate. We note, however, that polarization effects are crucial for divalent cations and only marginally important for monovalent cations (21). Therefore, although caution is warranted in accepting absolute values of binding free energies, qualitative trends and structural properties are likely to be well described by our approach.
Our results suggest that the protonation state is different in the sodium-and potassium-bound states. We show that although the protonation state controls potassium ion selectivity and stoichiometry, selectivity for sodium appears to be determined by the steric constraints imposed by the binding site geometry and independent from protonation state. Our simulations also suggest the presence of a previously unknown cytoplasmic binding site for anions that helps control the change in protonation state of Asp 926 , a crucial residue in forming the sodium-specific ion-binding site III. Furthermore, we show that access to this cytoplasmic anion-binding site is largely controlled by conformational changes in the cytoplasmic loops between transmembrane (TM) 4 helix 6 (TM6)-TM7 and TM8-TM9, which are different in sodium-bound E1P versus potassium-bound E2P states.

Results and discussion
The side chain carboxylate oxygens of five acidic amino acids (Glu 327 , Glu 779 , Asp 804 , Asp 808 , and Asp 926 ) are critical in coordinating Na ϩ and K ϩ in the Na ϩ /K ϩ -ATPase. Asp 926 is located on the TM8, and Glu 327 is located in TM4 (12), whereas the other acidic residues are located on TM5 and TM6 and are considerably closer to each other (Fig. 1). These five acidic residues form three binding sites: sites I and II are formed by Glu 327 , Glu 779 , Asp 804 , and Asp 808 along with backbone carbonyl moieties of TM4 residues and can be occupied by either sodium or potassium ions, and site III is formed by Asp 926 , Asp 808 , and oxygen atoms of TM5 residues and is only occupied in the E1P state by Na ϩ (12,13). As in previous computational studies of the Na ϩ /K ϩ -ATPase (19), we first tried to predict the protonation state of titratable residues using the PROPKA algorithm, an empirical method that relates desolvation effects and intraprotein interactions to the positions and chemical nature of residues proximate to the titratable sites (22). According to these calculations, only Asp 804 and Glu 779 are protonated in E1P, but for the E2P state, all are protonated except Asp 808 (supplemental Fig. S1). However, our molecular dynamics simulations initialized from these configurations produced, after 100 ns of equilibration, configurations inconsistent with these predictions: self-consistency tests applying PROPKA to the snapshots extracted from these equilibrated trajectories predicted protonation states significantly different from those in the initial structure (supplemental Fig. S1). Unsatisfied by these inconsistent results, we turned to more extensive, but more accurate, simulation-based methods to determine the protonation states for both sodiumbound E1P and potassium-bound E2P states.

Site III has a high affinity for both sodium and potassium ions when Asp 926 is deprotonated
With five acidic binding residues, systematic comparison of the 2 5 ϭ 32 possible protonation combinations for both the sodium-bound E1P and potassium-bound E2P states of Na ϩ / K ϩ -ATPase is a challenging task, requiring molecular dynamics simulations for a quarter of a million atoms for each possible protonation state. Fortunately, several protonation states could be eliminated. Our initial simulations of 210 ns, in which none of the acidic binding residues were protonated, displayed significant ion binding conformational inconsistency from crystal structures. In these simulations, both sodium and potassium ions relocated preferentially to site III (Fig. 1, C and D), whereas all crystal structures of potassium-bound states show that only sites I and II are occupied by potassium ions (14 -17). Because the carboxylate moiety of Asp 926 interacts exclusively with the cation bound to site III, we hypothesized that the protonation of this residue in the potassium-bound E2P state might result in detachment of the potassium ion from site III. To test this hypothesis, we continued the simulation after protonation of Asp 926 . Within several nanoseconds, the potassium ion bound in site III migrated toward sites I and II. We furthermore used metadynamics simulations (see "Experimental procedures") to calculate potential of mean force (PMF) along two collective variables (23) characterizing the bound ion density; the results revealed that protonation of Asp 926 leads to a more favorable binding of K ϩ to binding sites I and II and not to site III (Fig. 1E). Occupation of site III by Na ϩ ions, in contrast, is consistent with the presence of such a cation in the crystal structure, thus leading us to conclude that Asp 926 should be deprotonated in the Na ϩ -bound E1P state but not in the K ϩ -bound E2P state. A similar conclusion was reached by Nissen and co-workers (5).

The sodium-bound E1P and potassium-bound E2P conformations favor distinct protonation states
Next, we evaluated the remaining 2 4 ϭ 16 protonation state possibilities for the acidic binding residues: Glu 327 , Glu 779 , Asp 804 , and Asp 808 . As a first step, we tested the possibility that all four residues are deprotonated or protonated in both the E1P and E2P conformations. Not surprisingly, in equilibration simulations of each case, the large charge density in the binding site resulted in the binding of counterions from the solution buffer to neutralize the binding site charge, leading us to reject both possibilities (data not shown). For the remaining 14 combinations (supplemental Table S1), we were able to filter the results of equilibrium MD simulations using structural criteria, namely the root mean square deviation of binding site ions and residues from crystal structures, the degree of fluctuation of ions in the binding site, and ion coordination numbers (supplemental Figs. S2 and S3). We stress that rejection of particular protonation states from further consideration is not arbitrary; it is objectively based on the root mean square deviation from the experimentally determined structure. Our comprehensive analysis resulted in the elimination of seven additional combinations in the sodium-bound E1P state and eight additional protonation combinations in the potassium-bound E2P state (supplemental Tables S2 and S3).
For the remaining six and seven protonation combinations in the potassium-and sodium-bound conformations, respectively, we used free-energy perturbation (FEP) methods to determine the relative free energy of each protonation state (see "Experimental procedures"). To perform these calculations efficiently, alchemical intermediates were constructed along a network of pathways bridging each protonation state in the potassium-and sodium-bound conformations (Fig. 2). Forward and backward transformations, each starting from independent equilibrated structures to avoid bias, were used to validate the accuracy of these calculations. In total, this amounted to 10 FEP calculations for the potassium-bound E2P state (five forward and five backward) and 12 FEP calculations for the sodiumbound E1P state (six forward and six backward). Moreover, in simulations for which the total charge of the system changes during FEP, a reference state is used to compensate for anisotropic electrostatic effects (see "Experimental procedures"). Predicted free-energy differences, ⌬G, for all FEP simulations are summarized in Fig. 2. The results show that only two residues of four (Asp 804 , Asp 808 , Glu 327 , and Glu 779 ) are likely to be protonated in the sodium-bound and potassium-bound conformations, although the preferred protonated residues are different in each case. Specifically, for the sodium-bound E1P state, residues Glu 779 and Asp 804 are most likely to be protonated, whereas for the potassium-bound E2P state, residues Glu 327 and Asp 808 are most likely to be protonated (Fig. 2). and P-domain, green. The ␤-subunit is shown in violet, and the ␥-subunit is shown in orange. B, a superposition of the binding site for sodium-bound E1P (yellow) and potassium-bound E2P (green). C and D, density of ions sampled in 200-ns simulations (blue mesh; all ionizable binding site residues are unprotonated) superimposed on crystal structures for potassium (C) and sodium (D) where the green and yellow spheres represent the position of potassium and sodium ions, respectively, found in published Na ϩ /K ϩ -ATPase structures. E, a two-dimensional PMF calculated from metadynamics simulations for one potassium ion in the E2P binding site when Asp 926 is protonated. Approximate location of sites I, II, and III are highlighted with circles. F, scheme representing the collective variables used to calculate the PMF (see "Experimental procedures"). Black spheres show the reference atoms (␣-carbons) used for the metadynamics simulation.

FEP-determined protonation states correctly reproduce experimentally measured ion binding affinities
We performed further alchemical FEP calculations (see "Experimental procedures") to determine the relative affinities of sodium and potassium ions for the binding site in E1P and E2P conformations for the most likely protonation states described in the previous section (Fig. 2). In these calculations, sodium ions bound in the E1P state are alchemically transformed into potassium ions; likewise, potassium ions bound in the E2P state are transformed into sodium ions (see "Experimental procedures"). The results correctly predict the ion-binding selectivity of E1P and E2P states, estimating a preference of ϳ6 kcal/mol for sodium over potassium ions in the E1P conformation and a preference of ϳ1 kcal/mol for potassium over sodium ions in the E2P conformation. The experimentally measured selectivities are 1.5 kcal/mol preference for sodium in E1P and 4.5 kcal/mol preference for potassium ions in E2P (3,24,25). Although these experimental measurements are apparent affinities and cannot be directly compared with our calculated affinities, nevertheless, it is noteworthy that the calculation of ion affinities by FEP qualitatively reproduces ion selectivity in both E1P and E2P states.

Molecular mechanism of selectivity in Na ؉ /K ؉ -ATPase
Although the results above provide a description of distinct differences in protonation states and ion-binding selectivity for the E1P versus E2P conformations, a more important question Figure 2. FEP calculations of relative free energies of binding site protonation states for the sodium-bound E1P (A) and potassium-bound E2P (B) conformations of the Na ϩ /K ϩ -ATPase are shown. Protonated residues are labeled for each state and indicated by a letter "H" in each schematic. n.d., not determined. The inset in B is a reminder that Asp 926 is protonated in all these protonation states (see text). One, Two, Three, or Four Protons denote that one, two, three, or four acidic residues in the binding site are protonated, respectively. C and D, blue circles show the relative free energy of each protonation state obtained from the corresponding top panels with error bars (see "Experimental Procedures" for error bar calculation approach). Red stars show, for each protonation state, the change in free energy incurred about transforming three sodium ions into three potassium ions (C) and two potassium ions to two sodium ions (D). remains: what are the structural mechanisms underlying these distinct differences between the electrostatic environments of E1P and E2P? To provide insights, we used alchemical freeenergy calculations to determine the affinity of Na ϩ versus K ϩ in the E1P and E2P binding sites as well as electrostatic decomposition analysis to discover the origins of pK a shift for the binding site residues in E1P and E2P states.

Selectivity for sodium ions in the E1P state arises from binding site volume constraints
For all possible protonation states that show low deviations from crystal structures, i.e. the seven protonation states for E1P state and six protonation states in E2P state shown in Fig. 2, A and B, we calculated the binding free energy of sodium versus potassium to determine which ions are more preferred in each of these protonation states (Fig. 2, C and D). Our FEP results predict that in the E1P conformation the affinity for sodium ions relative to potassium ions does not depend on the specific protonation state. Regardless of protonation state, the E1P state always has higher affinity for sodium than potassium ions ( Fig.  2C and supplemental Fig. S4). Seemingly, only two reasons can justify this systematic increase in free energy upon change of sodium into potassium: a lack of ligands (water or protein oxygens) leading to an incomplete shell of solvation for the potassium ion or unfavorable repulsive van der Waals interactions between the binding pocket and the bigger (compared with sodium) solvated ion shell of potassium.
We found that the total coordination number increases as sodium ions are alchemically transformed to potassium ions, reaching values typical of a potassium ion in solution (supplemental Fig. S5); therefore, it is unlikely that potassium ions are suffering from lack of coordinating ligands. Conversely, the number of coordinating water molecules also increases as sodium ions are transformed into potassium ions (supplemental Fig. S5). These observations suggest that the E1P binding site cannot accommodate three larger solvated potassium ions. Accordingly, by analyzing separately each contribution to the average potential energy in two distinct simulations of the sodium-and potassium-bound states, we concluded that Lennard-Jones interactions have a destabilizing role (by as much as 8 kcal/mol) in passing from the first to the second case. We note that the same energy difference for an ion in aqueous environment has the opposite sign (it is approximately Ϫ2.5 kcal/mol; see supplemental Table S4 for further details). Although these enthalpy estimates are notoriously characterized by large statistical uncertainties (and despite the compensatory contribution expected from the entropic term), the Lennard-Jones energy difference is relatively large and suggestive of a significant role of steric interactions in ion selectivity. We thus surmise that the free-energy difference between the two states can be explained by the increased size of the solvated ion in the environment of a quasirigid pocket.
This result is consistent with the experimental finding that the E2P ion-binding site can accommodate acetamidinium, whereas the E1P ion-binding site does not, mostly likely due to the size and hydration structure of this organic cation (26). Furthermore, a study by Jorgensen et al. (27) shows that sodium ion-dependent enzyme phosphorylation is reduced to a greater degree with mutations that change side chain volume but preserve total charge of the side chain moiety, i.e. mutations E327D, E779D, and D804E, compared with mutations that change total negative charge but preserve side chain volume, i.e. mutations E327Q, E779Q, and D804N (supplemental Table S5). This suggests that the E1P state is sensitive to side chain volume changes, supporting a mechanism where binding site volume controls selectively.
Unlike the E1P state, the higher affinity of the E2P state for potassium ions compared with sodium ions appears to be more dependent on the specific protonation state. In this case, the looser steric restrictions allow either two sodium or potassium ions to bind; however, the clear trend is that, as pH decreases, the relative selectivity for potassium ions increases (Fig. 2D). This prediction is consistent with well established experimental observations (19,20). Hence, we surmise that the protonation state is a key factor in determining the selectivity of the E2P conformation. Next, we examine how the E2P state protonation state may be coupled to its conformational state.

A cytoplasmic binding site for anions controls the protonation state for Asp 926
To understand how protonation state is controlled in the Na ϩ /K ϩ -ATPase, we used electrostatic potential decomposition to investigate how the protein and/or its environment affects pK a of the five acidic binding site residues (see "Experimental procedures"). The goal of this postsimulation analysis is to qualitatively identify structural features that might contribute to differences in the electrostatic environment between E1P and E2P states and thereby give rise to pK a shifts of the acidic binding site residues.
Briefly, for the most probable protonation state in E1P and E2P, we calculated the electrostatic potential of placing a negative charge density on the region occupied by side chains of residues Asp 804 , Asp 808 , Asp 926 , Glu 327 , and Glu 779 (see "Experimental procedures"). In this analysis, relatively positive electrostatic potentials favor deprotonation, whereas relatively negative electrostatic potentials favor protonation. By comparing the electrostatic potentials between E1P and E2P, we can then determine whether conformational change tends to favor protonation or deprotonation. Although calculation of electrostatic potential does not provide a complete picture of the forces determining residue protonation, it nonetheless appears to be a reasonable approximation because the tendency for the five acidic residues to become protonated changes from the E2P to E1P conformation ( Fig. 3A; entire system) and is qualitatively similar to the predictions arising from the accurate FEP calculations of ⌬G (Fig. 2, A and B). For example, both sets of calculations predict that Glu 327 is more likely to become deprotonated, Glu 779 is more likely to become protonated, Asp 804 is more likely to become protonated, Asp 808 is more likely to become deprotonated, and Asp 926 is more likely to become deprotonated in the E2P state. The advantage of using electrostatic potential is the possibility to decompose the electrostatic potential into individual contributions from each component (protein, water, lipid, and ions).
Comparing the electrostatic potential generated by the entire system with the potential generated only by the protein reveals that the protein follows the same protonation/deprotonation patterns as those predicted by the entire system for all binding site residues in both E1P and E2P states except Asp 926 in the E2P state (Fig. 3A, blue box). The values for protein are, in general, more negative than the entire system because the total charge of the protein is highly negative and is biasing the potential toward negative values and, therefore, favoring residue protonation. However, only for Asp 926 in the E2P state does the protein favor deprotonation of this residue more than the system as a whole (see "Experimental procedures"). Therefore, protein conformational change is affecting electrostatic environment of Asp 926 differently than for the other acidic residues in the ion-binding pocket.
Asp 926 is of particular importance in the function of Na ϩ / K ϩ -ATPase; as we showed earlier, Asp 926 needs to be protonated in the E2P state to facilitate proper potassium binding to sites I and II and deprotonated in the E1P to allow sodium ion to bind site III. Although the protein conformational changes in E1P and E2P states seem to be responsible for protonation/ deprotonation of most binding site residues, environmental . Identification of a putative cytoplasmic anion-binding site for Na ؉ /K ؉ -ATPase. A, average electrostatic potential due to the protonation of key residues in the Na ϩ /K ϩ -binding site for both the E1P and E2P states. For each binding site residue, the bars show the electrostatic potential calculated using the entire system (yellow) and the subset of protein atoms (red) without considering contributions from that specific residue (see "Experimental procedures"). Error bars denote standard deviations. B, comparison of the relative contributions of lipids, water, solvent cations (K ϩ ), and anions (Cl Ϫ ) to the electrostatic potential of key binding residues in the E2P conformation. In each case, the maximum electrostatic potential values are subtracted to properly compare the relative contributions of each component (see supplemental Table S12 for values before substraction). C, molecular representation of the cytoplasmic anion-binding site in Na ϩ /K ϩ -ATPase. Only the ␣-subunit is shown for clarity. The bound chloride ion is shown as a magenta sphere. D, the distance between the TM6-TM7 loop and the TM8-TM9 loop (top panel) and chloride ion density near the cytoplasmic anion-binding site (lower panel) are compared between the E1P and E2P states. The distance is calculated using the C ␦ atom of Glu 821 and backbone nitrogen of Arg 933 . The chloride density is calculated by counting the number of chloride ions within 7 Å of the backbone nitrogen of Arg 933 . E, a cavity (blue-gray) extends from Asp 926 (shown in "green" licorice) to the cytoplasm in the E2P conformation (right) but is blocked in the E1P state (left). The program fpocket (51) was used for cavity calculations (see "Experimental procedures"). F, schematic representation of a proposed mechanism of TM6-TM7 lid opening, leading to the binding of a chloride ion, protonation of Asp 926 , and protonation rearrangement of other binding residues. factors seem to be controlling the protonation state of Asp 926 . Indeed, by decomposing the electrostatic potential into contributions arising from protein, water, lipid, cations (K ϩ for E2P), and anions (Cl Ϫ ), we found that the electrostatic interactions favoring protonation of Asp 926 in E2P state are largely due to the influence of anions ( Fig. 3B and supplemental Tables S11 and S12).
To understand how anions influence the pK a of Asp 926 , the anion density was compared between the E1P and E2P conformations. Calculation of a density difference map near this ionbinding residue (supplemental Fig. S7) reveals a putative cytoplasmic binding site for anions composed of the side chains of residues Asn 935 and Gln 940 and three backbone NH groups from residues Arg 933 , Arg 934 , and Asn 935 in the loop between TM8 and TM9 (Fig. 3C). In our simulations of the E2P conformation, these residues coordinate a chloride ion on the cytoplasmic face of the enzyme directly below Asp 926 (Fig. 3, C and  E). Intriguingly, we also found that this binding site has a structure similar to that of crystallographically determined chloridebinding sites in other proteins (28, 29) (supplemental Fig. S8). Moreover, our simulations show that in the E1P conformation access to the anion-binding site is blocked by salt bridges formed among arginine, aspartate, and glutamate residues in the TM6-TM7 loop and residues in the TM8-TM9 loop in conjunction with Glu 1011 , Tyr 1015 , and Tyr 1016 from the C terminus (supplemental Fig. S9). These observations suggest that the TM6-TM7 loop acts as a lid that gates the chloride-binding site, allowing access to anions in the E2P conformation and obscuring access in the E1P conformation, thus controlling the protonation state of Asp 926 and potentially determining the selectivity of the Na ϩ /K ϩ -ATPase.
Similar to the E1P state, the crystal structure of the potassiumbound E2P state shows close proximity of the TM6-TM7 loop and the TM8-TM9 loop (supplemental Fig. S10) with Glu 821 from the TM6-TM7 loop interacting with backbone nitrogen atoms of residues Arg 933 , Arg 934 , and Asn 935 ; hence, the crystal structures do not reveal the anion-binding site. However, simulations predict dissociation of this flexible TM6-TM7 loop (Fig. 3C) and association of the TM8-TM9 loop with a chloride ion soon thereafter (supplemental Movie 1) only in the E2P state (Fig. 3D). We also noticed that, after disassociation of TM6-TM7 loop, site III becomes hydrated via a cytoplasmic cavity that extends from the cytoplasm up to residue Asp 926 (Fig. 3E). In contrast, for the E1P state, this cytoplasmic cavity is truncated; hence, Asp 926 is not water-accessible from the cytoplasm. To confirm that the TM6-TM7 loop disassociation and anion binding are not a singular event and are statistically significant, we initiated five independent simulations of the potassium-bound E2P state started from the crystal structure conformation and continued for over 200 ns each (supplemental Fig. S11). The results show that TM6-TM7 loop always dissociates from TM8-TM9 loop, although it is possible that they reassociate. However, the overall dissociated-to-associated ratio is 4.6:1 (supplemental Fig. S11). Furthermore, simulations show that the chloride ion density increases near the cytoplasmic anion-binding site (supplemental Fig. S11) whenever the TM6-TM7 loop is disassociated from the TM8-TM9 loop.
To directly compare the sodium-bound E1P state with the potassium-bound E2P state, we were limited to the use of crystal structures obtained from Sus scrofa because the only crystal structures available for the sodium-bound E1P state are from this organism (12,13). Nevertheless, to further investigate the presence of the anion-binding site for the E2P state, we performed additional simulations of the high-resolution (2.4-Å) crystal structure of K ϩ -bound E2P (PDB code 2ZXE) from Squalus acanthias using our FEP-predicted protonation state for the S. scrofa E2P conformation (residues Glu 334 , Asp 815 , and Asp 933 in the 2ZXE crystal structure are protonated, which correspond to residues Glu 327 , Asp 808 , and Asp 926 in the 3KDP crystal structure, respectively). Similar to the low-resolution (3.5-Å) structure (PDB code 3KDP), the high-resolution crystal structure of K ϩ -bound E2P (PDB code 2ZXE) shows a glutamate residue side chain (Glu 828 in the 2ZXE structure, which corresponds to Glu 821 in the 3KDP structure) from the TM6-TM7 loop occupying the anion site (supplemental Fig. S12). Simulations started from this conformation also showed dissociation of this flexible loop (supplemental Fig. S13), although the dissociation time (ϳ400 ns) is slower compared with the low-resolution 3KDP structure. Furthermore, the chloride density increases soon after the TM6-TM7 loop is dissociated from the TM8-TM9 loop (supplemental Fig. S13), supporting the existence of the cytoplasmic anion-binding site.
For completeness, we also compared simulated ion-binding poses for a low-resolution (4.3-Å) structure of the E1P state (PDB code 4HQJ) and a high-resolution (2.8-Å) structure (PDB code 3WGU) (see supplemental Fig. S14). Similar coordination arrangements were found in both sets of simulations.
Poulsen et al. (5) have previously suggested that a cavity among TM5, TM7, and TM8 provides cytoplasmic access in the E2P conformation for protonation of Asp 926 and proton leak currents (30,31). Cytoplasmic access to this cavity was also postulated to be blocked in the E1P conformation by the C-terminal hydrogen-bond network (5). However, their simulations did not provide insight into why accessibility of this cavity could change the pK a of Asp 926 . Although not negating the importance of hydrogen bonding in controlling cytoplasmic accessibility of this cavity, our simulations suggest that anion site occupancy is the major factor determining the change in local electrostatic environment that controls the Asp 926 pK a change between the E1P and E2P conformations.

Inherited mutations and Hofmeister effects support the hypothesis of a cytoplasmic anion-binding site
The anion-binding event observed in our simulations and its energetic characterization provide a viable explanation for the conformation-dependent pattern of protonated residues. As the system relaxes from the E1P state to the E2P state, the TM6-TM7 lid opens up, and a chloride ion binds to the cytoplasmic anion-binding site (Fig. 3F). The consequences of this binding are 2-fold. (i) The electrostatic potential in the region of space occupied by Asp 926 is modified in such a way as to favor the protonated form of this residue, and (ii) the motion of the lid results in the formation of an aqueous cavity joining Asp 926 with the cytoplasm (Fig. 3F).
One of the most immediate predictions of this model is that the structural integrity and/or the sequence identity of the TM6-TM7 and TM8-TM9 loops should be required for proper function of the pump. Intriguingly, inherited mutations of the ATPA2 gene associated with hemiplegic migraine have been mapped to the TM8-TM9 loop region (7,32). Indeed, a singleresidue polymorphism, R933P, identified in these studies has been shown to dramatically disable the activity of the Na ϩ /K ϩ -ATPase (5,34,35). The mutation R933A has significant but less severe effects on activity (36), consistent with the loss of the backbone NH in R933P but not R933A. To further probe the effects of this mutation on anion binding, we performed an ϳ70-ns simulation of E2P R933P starting from a chloridebound conformation. The chloride ion leaves the binding site within 3 ns and remains dissociated (supplemental Movie 2). These observations do not exclude other explanations for the functional importance of TM8-TM9 loop residues, in particular the prevailing model that posits that TM6-TM7 loop mutations and C-terminal tyrosine mutations disrupt a hydrogenbonding network critical to the C-terminal cytoplasmic ion pathway (5). Our findings are remarkably consistent with this mode, but with the slight difference in interpretation that C-terminal tyrosine interactions with the Arg 933 side chain are crucial in anchoring the putative cytoplasmic anion-binding site.
Additional data consistent with the involvement of anions in the function of the Na ϩ /K ϩ -ATPase come from extensive studies characterizing a Hofmeister effect whereby chaotropic ions stabilize the E1P state (25,37) and slow the E1P-E2P interconversion rate (38). Although previous studies have ascribed the Hofmeister effect in the Na ϩ /K ϩ -ATPase and the sarcoplasmic reticulum Ca 2ϩ -ATPase to modified electrostatic interactions at the membrane/enzyme/cytoplasm interface (38 -40), the existenceofananion-bindingsitenearthisinterfaceisnotinconsistent with these data. Future simulation work will help to test this hypothesis by predicting binding affinities for different anionic chemical species.
We also note that the presence of Glu 821 bound to this site in the crystal structure of the E2P state is not inconsistent with the existence of an anion-binding site. Crystallographic conformations are often biased toward low-energy and/or low-entropy configurations. Glu 821 is in a region that is disordered in the crystal and likely highly mobile in the physiological state.
While this manuscript was under review, a study by Rui et al. (41) was published showing that protonation states for E1P and E2P are different in the binding aspartate residues. They found that Asp 804 is protonated in the E1P state and that Asp 808 and Asp 926 are protonated in the E2P state in agreement with our findings. Unlike our study, the authors assumed a protonated state for all glutamate residues without performing FEP calculations. Our FEP calculations determined protonation states for glutamate residues and showed that Glu 779 is protonated in the E1P state and that Glu 327 is protonated in the E2P state. Furthermore, although Rui et al. (41) attributed ion selectivity to protonation differences, our findings suggest that the selectivity in the E2P state is particularly protonation-dependent, whereas binding site volume constraints control selectivity in the E1P state irrespective of protonation state.

Conclusion
The mechanism by which ion selectivity is achieved in Na ϩ / K ϩ -ATPase is a long-standing question. The structural similarity of the sodium/potassium ion-binding site in both sodiumand potassium-occluded conformations suggests that the pump utilizes subtle mechanisms for controlling selectivity. Here, we have conducted a large-scale simulation study of Na ϩ / K ϩ -ATPase in both the E1P and E2P conformations and determined that the protonation state is likely to be different for sodium-bound versus potassium-bound states. Our modeling suggests that the potassium-bound E2P conformation is particularly sensitive to the protonation state, whereas for the sodium-bound E1P conformation, van der Waals interactions at the binding sites are more important. The change in protonation state for Asp 926 is particularly crucial for sodium and potassium ion selectivity and binding stoichiometry. Close inspection of our simulation results reveals a putative, previously unknown cytoplasmic anion-binding site whose occupancy promotes Asp 926 protonation by changing local electrostatic potential and increasing cytoplasmic accessibility of a cavity in the protein that extends to Asp 926 . Access to the cytoplasmic anionbinding site is itself controlled by the TM6-TM7 loop whose disposition is different in E1P versus E2P states and likely coupled to larger domain motions in the catalytic cycle of the pump triggered by ATP binding and hydrolysis.
We emphasize that the model outlined above provides specific predictions testable by further joint experimental and simulation studies. One key prediction is that mutation or truncation of residues in the TM8-TM9 loop should significantly disturb the pump's activity via changes in orientation and/or affinity of the putative cytoplasmic anion-binding site. Another is that the effects of specific Hofmeister ions should correlate with affinities to the cytoplasmic anion-binding site. Future studies will be able to directly test these predictions as well as to probe in detail the requisite conformational motions involved in transitions between E1 and E2 states of Na ϩ /K ϩ -ATPase.

System preparation
Molecular structures for the sodium-and potassium-bound Na ϩ /K ϩ -ATPase were obtained from the PDB: codes 4HQJ (E1P-ADP state) and 3KDP (E2P state). These structures were chosen because they share the same sequence and are obtained from the same organism, ␣ 1 and ␤ 1 isoforms of pig kidney. The CHARMM-GUI web interface (42) was used to insert the protein into a 1-palmytoyl-2-oleoyl-sn-glycero-3-phosphatidylcholine lipid bilayer and surrounded by TIP3P water molecules, Cl Ϫ , and Na ϩ or K ϩ ions for the E1P or E2P state, respectively. We introduced a phosphorylated aspartic acid at position 369 by using parameters present in the CHARMM27 force field for aspartate and phosphate moieties. The total number of particles is 233,660 and 261,581 for the E1P and E2P systems, respectively. A summary of system specifications is shown in Table 1.

Molecular dynamics equilibration simulations
The CHARMM27 force field was used for protein and water; CHARMM36 was used for lipids. Periodic boundary conditions were used for all of the MD simulations, and the electrostatic potential was evaluated using the particle mesh Ewald method (43) with grid spacing of 1.2 Å and non-bonded cutoff of 11 Å. The lengths of all bonds containing hydrogen were constrained with the SHAKE/RATTLE algorithm (44). The system was maintained at a temperature of 300 K and pressure of 1 atm using a Langevin thermostat and barostat as implemented in the MD code NAMD2. 10 (45). The reversible reference system propagator algorithms (rRESPA) multiple time step method was used (46) Tables S2 and  S3.

Metadynamics simulations
Initially, the equilibrated structures in which none of the acidic residues in the binding site were protonated were used for metadynamics simulations to study binding site preferences of Na ϩ and K ϩ ions. Then, Asp 926 was protonated in the E2P state to compare the results with the deprotonated Asp 926 case (see Fig. 1 and the main text). In all cases, two distanceZ collective variables (23) were used to build the PMF and to confine one Na ϩ or K ϩ ion in a rectangular box of 30 ϫ 20 Å 2 inside the binding site. Potential hills with a height of 0.05 kcal/mol and a width of 0.25 Å were added every 200 steps for a total simulation time of 70 ns using grid-based storage of bias potentials.

FEP calculations to identify the lowest-energy protonation state
FEP calculations were performed by linearly interpolating the Hamiltonians of the protonated and deprotonated systems, (1 Ϫ )U prot ϩ U deprot , with the progress variable varying between 0 and 1. For each transformation, trajectories for different values of (windows), with each window simulated for 10 -60 ns, were collected. Additional harmonic potentials (with a force constant 50.0 kcal/mol/Å 2 applied to side chain atoms of the decoupling and coupling residues) were added to the potential energy function to keep the particles of the decoupling and coupling amino acids in close proximity. These additional interactions have negligible contribution to the free-energy estimates, as confirmed by superposition of final snapshots of FEP transformations with final snapshots of the equilibrium simulations to ensure that these constraints do not constrain the conformational space (data are not shown), while allowing contiguous windows to sample overlapping regions of the conformational space, thereby increasing the rates of convergence of the free-energy estimates. Subsequent analysis of trajectories was limited to the stable part of each window, i.e. the one with no detectable drift in the potential energy. Simulation times for all FEP runs are shown in supplemental Table S6. Initially, 20-window FEPs were used; each window followed from the equilibrated structure of the previous window to increase the rate of convergence. Then, in selected cases (when the ⌬⌬G was affected by large uncertainties), we used 40 FEP windows, all started with the same equilibrated structure obtained from at least 100 ns of unbiased equilibration MD, to estimate the free energy. In all cases, the results of 40-window FEPs were more accurate than 20-window FEPs, and reassuringly, they were always in qualitative agreement with those obtained with 20-window FEPs (supplemental Tables S7 and S8). Therefore, we performed 40-window FEPs for all transformations and reported the final results presented in Fig. 2 based on these FEPs. In both E1P and E2P states, all alchemical transformations except two entailed the transfer of a proton from one binding site residue to another. In these cases, the overall chemical composition of the system is unchanged during the transformation, and the estimated free-energy difference is directly comparable with the experiments. In the remaining cases, where the total charge of the system changes during FEP, an auxiliary thermodynamic cycle is considered in which the FEP of a proton dissociation reaction from a three-residue peptide is calculated in solution. Then, this FEP is subtracted from the FEP of protein to compensate charge-changing effects (see below for more details). For FEP reactions in which the total charge of the system is unchanged during the FEP alchemical transformation, a thermodynamic cycle was not used. Indeed, if a proton is transferred from one aspartate to another (or similarly from a glutamate to another), then the ⌬G of the reference reaction identically vanishes. Moreover, because aspartic and glutamic acids have similar pK a values in bulk water, the ⌬G of the reference reaction is negligible.

FEP calculations for reference (solution) reactions
The reference proton dissociation reactions for Glu and Asp amino acids in solution were characterized by using a tripeptide, Ile-Asp-Leu for Asp protonation and Pro-Glu-Ile for Glu protonation. These sequences were chosen on the basis of the neighboring residues of Asp 804 and Glu 779 in the Na ϩ /K ϩ -ATPase. Both systems were solvated with TIP3P water molecules and minimized before FEP runs. Similar to other FEP calculations, harmonic restraints were applied to keep the particles of decoupling and coupling amino acids in close proximity. 40 windows were used each with 10-ns equilibration time. All other parameters are the same as described above. See supplemental Table S9 for the FEP values.
Instead of using a reference reaction to cancel charge changing effects during FEP calculations, one may use protonation/ deprotonation of a side chain in the protein that is far from binding site; therefore, the charge of the system will stay the same during the FEP. We tested the accuracy of our reference Glu 327 -Glu 779 -Asp 808 transformation, the result is ambiguous because the error bar for this FEP reaction is more than the FEP value; therefore, the comparison of two approaches could not be achieved. Nevertheless, results of both FEP methods may be interpreted as that the Glu 327 -Glu 779 and Glu 327 -Glu 779 -Asp 808 protonation states have similar free energies, and both have higher energy than Glu 327 -Asp 808 . Therefore, these will not change the predicted lowest free-energy protonated state for the E2P conformation shown in Fig. 2 of the main text.

FEP calculations to determine Na ؉ versus K ؉ selectivity
To calculate the free-energy difference between the sodiumbound and potassium-bound states, the Lennard-Jones parameters of the bound ion were altered to produce a smooth transformation from Na ϩ to K ϩ (supplemental Fig. S6). 20 FEP windows were used, and each window was equilibrated for 5-10 ns with the free-energy difference between adjacent windows estimated according to the Zwanzig relation: ⌬G(A 3 B) ϭ Ϫk B T ln͗exp(Ϫ(E B Ϫ E A )/k B T)͘ A (47). To ensure its accuracy, we tested this protocol to calculate the difference in solvation free energy between sodium and potassium. The result (Ϫ19.8 Ϯ 0.1 kcal/mol) is within 1 kcal/mol of previously reported values (Ϫ20.7 Ϯ 0.1 kcal/mol) and closer to the experimental value (Ϫ17.5 kcal/mol (48)). These FEP values are multiplied by 2 to account for changing two potassium ions to two sodium ions in the E2P state and multiplied by 3 to account for changing three sodium ions to three potassium ions in the E1P state. For these alchemical transformations, for following thermodynamic cycle is used.

Analysis of FEP simulations
The multistate Bennett acceptance ratio method was used for free-energy estimation as implemented in the pymbar Python package (49,50). Because our FEP calculations to determine the preferred protonation state show good energy overlap not only between adjacent but also between more distant windows, we used the multistate Bennett acceptance ratio estimator, which integrates samples from all thermodynamic intermediates. This estimator is expected to produce superior estimates compared with exponential averaging and/or the (two-state) Bennett acceptance ratio, which only consider adjacent windows.

Electrostatic potential calculations
To gauge the relative propensity for the deprotonated state of residues Asp 926 , Glu 327 , Glu 779 , Asp 804 , and Asp 808 , the electrostatic potential energy of a negative charge density located in the regions occupied by each of these residues was estimated. This quantity can be interpreted as the electrostatic component of the free-energy difference on changing from the protonated to deprotonated state. We estimated this energy using the following equation.
where (r) is the density of the atoms belonging to side chains of residues Asp 926 , Glu 327 , Glu 779 , Asp 804 , or Asp 808 ; ͗V ext (r)͘ t is the electrostatic potential, averaged over the entire trajectory, generated by all the charges in the system except those contributing to (r). Because the system is periodic, the long-range part of the potential has been evaluated using the particle mesh Ewald method (43) with grid spacing of 1.0 Å and non-bonded cutoff of 11 Å. In this approach, in which the Poisson equation is solved in the reciprocal space, the k ϭ 0 term is customarily neglected as it is divergent. Therefore, ͗V ext (r)͘ t is defined only up to an additive constant; we use this property in supplemental Tables S11 and S12 to define a convenient baseline for the interaction energies. (r) is estimated from a three-dimensional histogram of atom positions using Gaussian kernels of width i equal to the van der Waals radius of each particle. In this equation, r indicates the position in space, r i (t) is the location of the ith side chain atom at time t, the sum runs over all the atoms of the binding site residue under consideration, and ͗ . . . ͘ t indicates the average over the MD trajectory. By virtue of the linearity of the Poisson equation, the external potential (and thus the electrostatic potential) can be decomposed in contributions arising from five distinct, non-overlapping sets of atoms: protein, water, lipid, solution (K ϩ or Na ϩ ) cations, and solution Cl Ϫ anions. The electrostatic potentials for each component and for all binding site residues are shown in supplemental Tables S11 and S12. Calculations in practice use grid-based storage of charge densities and external potentials.

Protein cavity detection
To identify pockets and pathways connecting the binding site to the intracellular milieu, we used the program fpocket (33,