Advertisement

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

  • Asghar M. Razavi
    Affiliations
    Department of Chemistry, Temple University, Philadelphia, Pennsylvania 19122
    Search for articles by this author
  • Lucie Delemotte
    Footnotes
    Affiliations
    Institute for Computational Molecular Science, Temple University, Philadelphia, Pennsylvania 19122

    Science for Life Laboratory, Department of Theoretical Physics, KTH Royal Institute of Technology, Stockholm 11428, Sweden
    Search for articles by this author
  • Joshua R. Berlin
    Affiliations
    Department of Pharmacology, Physiology and Neuroscience, New Jersey Medical School, Rutgers University, Newark, New Jersey 07103
    Search for articles by this author
  • Vincenzo Carnevale
    Correspondence
    To whom correspondence may be addressed: Institute for Computational Molecular Science, Temple University, SERC 704C, N. 12th St., Philadelphia, PA 19122. Tel.: 215-204-4214;.
    Affiliations
    Department of Chemistry, Temple University, Philadelphia, Pennsylvania 19122

    Institute for Computational Molecular Science, Temple University, Philadelphia, Pennsylvania 19122
    Search for articles by this author
  • Vincent A. Voelz
    Correspondence
    To whom correspondence may be addressed: Dept. of Chemistry, Temple University, Beury 240, 1901 N. 13th St., Philadelphia, PA 19122. Tel.: 215-204-1973;
    Affiliations
    Department of Chemistry, Temple University, Philadelphia, Pennsylvania 19122

    Institute for Computational Molecular Science, Temple University, Philadelphia, Pennsylvania 19122
    Search for articles by this author
  • Author Footnotes
    1 Recipient of funding from European Union Seventh Framework Program 979 “Voltsens” Grant PIOF-GA-2012-329534.
Open AccessPublished:June 06, 2017DOI:https://doi.org/10.1074/jbc.M117.779090
      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.

      Introduction

      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 (
      • Skou J.C.
      The identification of the sodium pump.
      ,
      • Skou J.C.
      The influence of some cations on an adenosine triphosphatase from peripheral nerves.
      ,
      • Jorgensen P.L.
      • Hakansson K.O.
      • Karlish S.J.
      Structure and mechanism of Na,K-ATPase: functional sites and their interactions.
      ). 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 (
      • Schwinger R.H.
      • Bundgaard H.
      • Müller-Ehmsen J.
      • Kjeldsen K.
      The Na, K-ATPase in the failing human heart.
      ,
      • Poulsen H.
      • Khandelia H.
      • Morth J.P.
      • Bublitz M.
      • Mouritsen O.G.
      • Egebjerg J.
      • Nissen P.
      Neurological disease mutations compromise a C-terminal ion pathway in the Na+/K+-ATPase.
      ,
      • Moseley A.E.
      • Williams M.T.
      • Schaefer T.L.
      • Bohanan C.S.
      • Neumann J.C.
      • Behbehani M.M.
      • Vorhees C.V.
      • Lingrel J.B.
      Deficiency in Na,K-ATPase α isoform genes alters spatial learning, motor activity, and anxiety in mice.
      ); thus, Na+/K+-ATPase is an important target for treatment of brain and heart conditions (
      • Morth J.P.
      • Poulsen H.
      • Toustrup-Jensen M.S.
      • Schack V.R.
      • Egebjerg J.
      • Andersen J.P.
      • Vilsen B.
      • Nissen P.
      The structure of the Na+, K+-ATPase and mapping of isoform differences and disease-related mutations.
      ,
      • Yatime L.
      • Buch-Pedersen M.J.
      • Musgaard M.
      • Morth J.P.
      • Lund Winther A.M.
      • Pedersen B.P.
      • Olesen C.
      • Andersen J.P.
      • Vilsen B.
      • Schiøtt B.
      • Palmgren M.G.
      • Møller J.V.
      • Nissen P.
      • Fedosova N.
      P-type ATPases as drug targets: tools for medicine and science.
      ).
      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 (
      • Post R.L.
      • Sen A.K.
      • Rosenthal A.S.
      A phosphorylated intermediate in adenosine triphosphate-dependent sodium and potassium transport across kidney membranes.
      ,
      • Albers R.
      Biochemical aspects of active transport.
      ,
      • Post R.L.
      • Kume S.
      • Tobin T.
      • Orcutt B.
      • Sen A.K.
      Flexibility of an active center in sodium-plus-potassium adenosine triphosphatase.
      ), here we focus on the sodium-bound E1P and potassium-bound E2P states (i.e. phosphorylated E1 and E2) for which crystal structures are available (
      • Kanai R.
      • Ogawa H.
      • Vilsen B.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of a Na+-bound Na+,K+-ATPase preceding the E1P state.
      ,
      • Nyblom M.
      • Poulsen H.
      • Gourdon P.
      • Reinhard L.
      • Andersson M.
      • Lindahl E.
      • Fedosova N.
      • Nissen P.
      Crystal structure of Na+, K+-ATPase in the Na+-bound state.
      ,
      • Laursen M.
      • Yatime L.
      • Nissen P.
      • Fedosova N.U.
      Crystal structure of the high-affinity Na+,K+-ATPase–ouabain complex with Mg2+ bound in the cation binding site.
      ,
      • Morth J.P.
      • Pedersen B.P.
      • Toustrup-Jensen M.S.
      • Sørensen T.L.
      • Petersen J.
      • Andersen J.P.
      • Vilsen B.
      • Nissen P.
      Crystal structure of the sodium-potassium pump.
      ,
      • Ogawa H.
      • Shinoda T.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of the sodium-potassium pump (Na+,K+-ATPase) with bound potassium and ouabain.
      ,
      • Shinoda T.
      • Ogawa H.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of the sodium-potassium pump at 2.4 Å resolution.
      ). One of the central features of ion transport by the Na+/K+-ATPase is a change in ion selectivity between E1 and E2 conformations (
      • Jorgensen P.L.
      • Hakansson K.O.
      • Karlish S.J.
      Structure and mechanism of Na,K-ATPase: functional sites and their interactions.
      ). The origin of this selectivity change has been a focus of Na+/K+-ATPase research since its discovery by Jens Christian Skou in 1957 (
      • Skou J.C.
      The influence of some cations on an adenosine triphosphatase from peripheral nerves.
      ). Extensive mutational studies have identified key residues involved in ion binding, which include five acidic residues, Asp804, Asp808, and Asp926, Glu327, and Glu779 (
      • Nyblom M.
      • Poulsen H.
      • Gourdon P.
      • Reinhard L.
      • Andersson M.
      • Lindahl E.
      • Fedosova N.
      • Nissen P.
      Crystal structure of Na+, K+-ATPase in the Na+-bound state.
      ), whose orientation, distance, and ion coordination have been proposed to be involved in determining selectivity (
      • Jorgensen P.L.
      • Hakansson K.O.
      • Karlish S.J.
      Structure and mechanism of Na,K-ATPase: functional sites and their interactions.
      ). The first crystal structure of the Na+/K+-ATPase, solved for the potassium-bound E2P state in 2007 (
      • Morth J.P.
      • Pedersen B.P.
      • Toustrup-Jensen M.S.
      • Sørensen T.L.
      • Petersen J.
      • Andersen J.P.
      • Vilsen B.
      • Nissen P.
      Crystal structure of the sodium-potassium pump.
      ), and later followed by several higher resolution structures (
      • Laursen M.
      • Yatime L.
      • Nissen P.
      • Fedosova N.U.
      Crystal structure of the high-affinity Na+,K+-ATPase–ouabain complex with Mg2+ bound in the cation binding site.
      ,
      • Ogawa H.
      • Shinoda T.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of the sodium-potassium pump (Na+,K+-ATPase) with bound potassium and ouabain.
      ,
      • Shinoda T.
      • Ogawa H.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of the sodium-potassium pump at 2.4 Å resolution.
      ) 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 Ca2+ pump, sarco/endoplasmic reticulum Ca2+-ATPase, a closely related P-type ATPase (
      • Møller J.
      • Nissen P.
      • Sørensen T.
      • le Maire M.
      Transport mechanism of the sarcoplasmic reticulum Ca2+-ATPase pump.
      ).
      Figure thumbnail gr1
      Figure 1Molecular structure of Na+/K+-ATPase. A, the α-subunit is shown in gray with its cytoplasmic domains colored as follows: A-domain, red; N-domain, blue; 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 Asp926 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.
      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. (
      • Yu H.
      • Ratheal I.M.
      • Artigas P.
      • Roux B.
      Protonation of key acidic residues is critical for the K+-selectivity of the Na/K pump.
      ) in which different protonation states of the potassium-bound E2P state were predicted to confer differences in ion selectivity. In this work, electrophysiological experiments confirmed computational predictions that potassium selectivity decreases with increasing pH, consistent with previous pH studies (
      • Milanick M.A.
      • Arnett K.L.
      Extracellular protons regulate the extracellular cation selectivity of the sodium pump.
      ). 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 sodium-bound E1P state were obtained in 2013 (
      • Kanai R.
      • Ogawa H.
      • Vilsen B.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of a Na+-bound Na+,K+-ATPase preceding the E1P state.
      ,
      • Nyblom M.
      • Poulsen H.
      • Gourdon P.
      • Reinhard L.
      • Andersson M.
      • Lindahl E.
      • Fedosova N.
      • Nissen P.
      Crystal structure of Na+, K+-ATPase in the Na+-bound state.
      ), 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 (
      • Li H.
      • Ngo V.
      • Da Silva M.C.
      • Salahub D.R.
      • Callahan K.
      • Roux B.
      • Noskov S.Y.
      Representation of ion–protein interactions using the Drude polarizable force-field.
      ). 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 Asp926, 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)
      The abbreviations used are: TM
      transmembrane
      MD
      molecular dynamics
      FEP
      free-energy perturbation
      PMF
      potential of mean force
      PDB
      Protein Data Bank.
      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 (Glu327, Glu779, Asp804, Asp808, and Asp926) are critical in coordinating Na+ and K+ in the Na+/K+-ATPase. Asp926 is located on the TM8, and Glu327 is located in TM4 (
      • Kanai R.
      • Ogawa H.
      • Vilsen B.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of a Na+-bound Na+,K+-ATPase preceding the E1P state.
      ), 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 Glu327, Glu779, Asp804, and Asp808 along with backbone carbonyl moieties of TM4 residues and can be occupied by either sodium or potassium ions, and site III is formed by Asp926, Asp808, and oxygen atoms of TM5 residues and is only occupied in the E1P state by Na+ (
      • Kanai R.
      • Ogawa H.
      • Vilsen B.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of a Na+-bound Na+,K+-ATPase preceding the E1P state.
      ,
      • Nyblom M.
      • Poulsen H.
      • Gourdon P.
      • Reinhard L.
      • Andersson M.
      • Lindahl E.
      • Fedosova N.
      • Nissen P.
      Crystal structure of Na+, K+-ATPase in the Na+-bound state.
      ). As in previous computational studies of the Na+/K+-ATPase (
      • Yu H.
      • Ratheal I.M.
      • Artigas P.
      • Roux B.
      Protonation of key acidic residues is critical for the K+-selectivity of the Na/K pump.
      ), 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 (
      • Li H.
      • Robertson A.D.
      • Jensen J.H.
      Very fast empirical prediction and rationalization of protein pKa values.
      ). According to these calculations, only Asp804 and Glu779 are protonated in E1P, but for the E2P state, all are protonated except Asp808 (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 sodium-bound E1P and potassium-bound E2P states.

      Site III has a high affinity for both sodium and potassium ions when Asp926 is deprotonated

      With five acidic binding residues, systematic comparison of the 25 = 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 (
      • Laursen M.
      • Yatime L.
      • Nissen P.
      • Fedosova N.U.
      Crystal structure of the high-affinity Na+,K+-ATPase–ouabain complex with Mg2+ bound in the cation binding site.
      ,
      • Morth J.P.
      • Pedersen B.P.
      • Toustrup-Jensen M.S.
      • Sørensen T.L.
      • Petersen J.
      • Andersen J.P.
      • Vilsen B.
      • Nissen P.
      Crystal structure of the sodium-potassium pump.
      ,
      • Ogawa H.
      • Shinoda T.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of the sodium-potassium pump (Na+,K+-ATPase) with bound potassium and ouabain.
      ,
      • Shinoda T.
      • Ogawa H.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of the sodium-potassium pump at 2.4 Å resolution.
      ). Because the carboxylate moiety of Asp926 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 Asp926. 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 (
      • Fiorin G.
      • Klein M.L.
      • Hénin J.
      Using collective variables to drive molecular dynamics simulations.
      ) characterizing the bound ion density; the results revealed that protonation of Asp926 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 Asp926 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 (
      • Poulsen H.
      • Khandelia H.
      • Morth J.P.
      • Bublitz M.
      • Mouritsen O.G.
      • Egebjerg J.
      • Nissen P.
      Neurological disease mutations compromise a C-terminal ion pathway in the Na+/K+-ATPase.
      ).

      The sodium-bound E1P and potassium-bound E2P conformations favor distinct protonation states

      Next, we evaluated the remaining 24 = 16 protonation state possibilities for the acidic binding residues: Glu327, Glu779, Asp804, and Asp808. 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 sodium-bound 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 (Asp804, Asp808, Glu327, and Glu779) 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 Glu779 and Asp804 are most likely to be protonated, whereas for the potassium-bound E2P state, residues Glu327 and Asp808 are most likely to be protonated (Fig. 2).
      Figure thumbnail gr2
      Figure 2FEP 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 Asp926 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).

      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 (
      • Jorgensen P.L.
      • Hakansson K.O.
      • Karlish S.J.
      Structure and mechanism of Na,K-ATPase: functional sites and their interactions.
      ,
      • Pedersen P.A.
      • Nielsen J.M.
      • Rasmussen J.H.
      • Jorgensen P.L.
      Contribution to Tl+, K+, and Na+ binding of Asn776, Ser775, Thr774, Thr772, and Tyr771 in cytoplasmic part of fifth transmembrane segment in R-subunit of renal Na,K-ATPase.
      ,
      • Klodos I.
      • Post R.L.
      • Forbush 3rd, B.
      Kinetic heterogeneity of phosphoenzyme of Na,K-ATPase modeled by unmixed lipid phases.
      ). 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 remains: what are the structural mechanisms underlying these distinct differences between the electrostatic environments of E1P and E2P? To provide insights, we used alchemical free-energy 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 pKa 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 (
      • Mahmmoud Y.A.
      • Kopec W.
      • Khandelia H.
      K+ congeners that do not compromise Na+ activation of the Na+, K+-ATPase: hydration of the ion binding cavity likely controls ion selectivity.
      ). Furthermore, a study by Jorgensen et al. (
      • Jorgensen P.L.
      • Rasmussen J.H.
      • Nielsen J.M.
      • Pedersen P.A.
      Transport-linked conformational changes in Na,K-ATPase. Structure-function relationships of ligand binding and E1-E2 conformational transitions.
      ) 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 (
      • Yu H.
      • Ratheal I.M.
      • Artigas P.
      • Roux B.
      Protonation of key acidic residues is critical for the K+-selectivity of the Na/K pump.
      ,
      • Milanick M.A.
      • Arnett K.L.
      Extracellular protons regulate the extracellular cation selectivity of the sodium pump.
      ). 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 Asp926

      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 pKa 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 pKa 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 Asp804, Asp808, Asp926, Glu327, and Glu779 (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 Glu327 is more likely to become deprotonated, Glu779 is more likely to become protonated, Asp804 is more likely to become protonated, Asp808 is more likely to become deprotonated, and Asp926 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).
      Figure thumbnail gr3
      Figure 3Identification 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 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 Glu821 and backbone nitrogen of Arg933. The chloride density is calculated by counting the number of chloride ions within 7 Å of the backbone nitrogen of Arg933. E, a cavity (blue-gray) extends from Asp926 (shown in “green” licorice) to the cytoplasm in the E2P conformation (right) but is blocked in the E1P state (left). The program fpocket (
      • Le Guilloux V.
      • Schmidtke P.
      • Tuffery P.
      Fpocket: an open source platform for ligand pocket detection.
      ) 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 Asp926, and protonation rearrangement of other binding residues.
      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 Asp926 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 Asp926 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 Asp926 differently than for the other acidic residues in the ion-binding pocket.
      Asp926 is of particular importance in the function of Na+/K+-ATPase; as we showed earlier, Asp926 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 factors seem to be controlling the protonation state of Asp926. 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 Asp926 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 pKa of Asp926, the anion density was compared between the E1P and E2P conformations. Calculation of a density difference map near this ion-binding residue (supplemental Fig. S7) reveals a putative cytoplasmic binding site for anions composed of the side chains of residues Asn935 and Gln940 and three backbone NH groups from residues Arg933, Arg934, and Asn935 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 Asp926 (Fig. 3, C and E). Intriguingly, we also found that this binding site has a structure similar to that of crystallographically determined chloride-binding sites in other proteins (
      • Lim H.H.
      • Shane T.
      • Miller C.
      Intracellular proton access in a Cl/H+ antiporter.
      ,
      • Ogawa H.
      • Qiu Y.
      • Philo J.S.
      • Arakawa T.
      • Ogata C.M.
      • Misono K.S.
      Reversibly bound chloride in the atrial natriuretic peptide receptor hormone-binding domain: possible allosteric regulation and a conserved structural motif for the chloride-binding site.
      ) (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 Glu1011, Tyr1015, and Tyr1016 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 Asp926 and potentially determining the selectivity of the Na+/K+-ATPase.
      Similar to the E1P state, the crystal structure of the potassium-bound E2P state shows close proximity of the TM6-TM7 loop and the TM8-TM9 loop (supplemental Fig. S10) with Glu821 from the TM6-TM7 loop interacting with backbone nitrogen atoms of residues Arg933, Arg934, and Asn935; 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 Asp926 (Fig. 3E). In contrast, for the E1P state, this cytoplasmic cavity is truncated; hence, Asp926 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 (
      • Kanai R.
      • Ogawa H.
      • Vilsen B.
      • Cornelius F.
      • Toyoshima C.
      Crystal structure of a Na+-bound Na+,K+-ATPase preceding the E1P state.
      ,
      • Nyblom M.
      • Poulsen H.
      • Gourdon P.
      • Reinhard L.
      • Andersson M.
      • Lindahl E.
      • Fedosova N.
      • Nissen P.
      Crystal structure of Na+, K+-ATPase in the Na+-bound state.
      ). 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 Glu334, Asp815, and Asp933 in the 2ZXE crystal structure are protonated, which correspond to residues Glu327, Asp808, and Asp926 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 (Glu828 in the 2ZXE structure, which corresponds to Glu821 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. (
      • Poulsen H.
      • Khandelia H.
      • Morth J.P.
      • Bublitz M.
      • Mouritsen O.G.
      • Egebjerg J.
      • Nissen P.
      Neurological disease mutations compromise a C-terminal ion pathway in the Na+/K+-ATPase.
      ) have previously suggested that a cavity among TM5, TM7, and TM8 provides cytoplasmic access in the E2P conformation for protonation of Asp926 and proton leak currents (
      • Li C.
      • Geering K.
      • Horisberger J.
      The third sodium binding site of Na,K-ATPase is functionally linked to acidic pH-activated inward current.
      ,
      • Yaragatupalli S.
      • Olivera J.F.
      • Gatto C.
      • Artigas P.
      Altered Na+ transport after an intracellular α-subunit deletion reveals strict external sequential release of Na+ from the Na/K pump.
      ). Cytoplasmic access to this cavity was also postulated to be blocked in the E1P conformation by the C-terminal hydrogen-bond network (
      • Poulsen H.
      • Khandelia H.
      • Morth J.P.
      • Bublitz M.
      • Mouritsen O.G.
      • Egebjerg J.
      • Nissen P.
      Neurological disease mutations compromise a C-terminal ion pathway in the Na+/K+-ATPase.
      ). However, their simulations did not provide insight into why accessibility of this cavity could change the pKa of Asp926. 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 Asp926 pKa 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 Asp926 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 Asp926 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 (
      • Morth J.P.
      • Poulsen H.
      • Toustrup-Jensen M.S.
      • Schack V.R.
      • Egebjerg J.
      • Andersen J.P.
      • Vilsen B.
      • Nissen P.
      The structure of the Na+, K+-ATPase and mapping of isoform differences and disease-related mutations.
      ,
      • Riant F.
      • De Fusco M.
      • Aridon P.
      • Ducros A.
      • Ploton C.
      • Marchelli F.
      • Maciazek J.
      • Bousser M.G.
      • Casari G.
      • Tournier-Lasserve E.
      ATP1A2 mutations in 11 families with familial hemiplegic migraine.
      ). Indeed, a single-residue polymorphism, R933P, identified in these studies has been shown to dramatically disable the activity of the Na+/K+-ATPase (
      • Poulsen H.
      • Khandelia H.
      • Morth J.P.
      • Bublitz M.
      • Mouritsen O.G.
      • Egebjerg J.
      • Nissen P.
      Neurological disease mutations compromise a C-terminal ion pathway in the Na+/K+-ATPase.
      ,
      • Tavraz N.N.
      • Friedrich T.
      • Dürr K.L.
      • Koenderink J.B.
      • Bamberg E.
      • Freilinger T.
      • Dichgans M.
      Diverse functional consequences of mutations in the Na+/K+-ATPase α2-subunit causing familial hemiplegic migraine type 2.
      ,
      • Weigand K.M.
      • Swarts H.G.
      • Russel F.G.
      • Koenderink J.B.
      Biochemical characterization of sporadic/familial hemiplegic migraine mutations.
      ). The mutation R933A has significant but less severe effects on activity (
      • Toustrup-Jensen M.S.
      • Holm R.
      • Einholm A.P.
      • Schack V.R.
      • Morth J.P.
      • Nissen P.
      • Andersen J.P.
      • Vilsen B.
      The C terminus of Na+,K+-ATPase controls Na+ affinity on both sides of the membrane through Arg935.
      ), 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 chloride-bound 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 hydrogen-bonding network critical to the C-terminal cytoplasmic ion pathway (
      • Poulsen H.
      • Khandelia H.
      • Morth J.P.
      • Bublitz M.
      • Mouritsen O.G.
      • Egebjerg J.
      • Nissen P.
      Neurological disease mutations compromise a C-terminal ion pathway in the Na+/K+-ATPase.
      ). Our findings are remarkably consistent with this mode, but with the slight difference in interpretation that C-terminal tyrosine interactions with the Arg933 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 (
      • Klodos I.
      • Post R.L.
      • Forbush 3rd, B.
      Kinetic heterogeneity of phosphoenzyme of Na,K-ATPase modeled by unmixed lipid phases.
      ,
      • Suzuki K.
      • Post R.
      Equilibrium of phosphointermediates of sodium and potassium ion transport adenosine triphosphatase. Action of sodium ion and Hofmeister effect.
      ) and slow the E1P-E2P interconversion rate (
      • Ganea C.
      • Babes A.
      • Lüpfert C.
      • Grell E.
      • Fendler K.
      • Clarke R.J.
      Hofmeister effects of anions on the kinetics of partial reactions of the Na+,K+-ATPase.
      ). Although previous studies have ascribed the Hofmeister effect in the Na+/K+-ATPase and the sarcoplasmic reticulum Ca2+-ATPase to modified electrostatic interactions at the membrane/enzyme/cytoplasm interface (
      • Ganea C.
      • Babes A.
      • Lüpfert C.
      • Grell E.
      • Fendler K.
      • Clarke R.J.
      Hofmeister effects of anions on the kinetics of partial reactions of the Na+,K+-ATPase.
      ,
      • Tadini-Buoninsegni F.
      • Moncelli M.R.
      • Peruzzi N.
      • Ninham B.W.
      • Dei L.
      • Nostro P.L.
      Hofmeister effect of anions on calcium translocation by sarcoplasmic reticulum Ca2+-ATPase.
      ,
      • Clarke R.J.
      Dipole-potential-mediated effects on ion pump kinetics.
      ), the existence of an anion-binding site near this interface is not inconsistent 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 Glu821 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. Glu821 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. (
      • Rui H.
      • Artigas P.
      • Roux B.
      The selectivity of the Na+/K+-pump is controlled by binding site protonation and self-correcting occlusion.
      ) was published showing that protonation states for E1P and E2P are different in the binding aspartate residues. They found that Asp804 is protonated in the E1P state and that Asp808 and Asp926 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 Glu779 is protonated in the E1P state and that Glu327 is protonated in the E2P state. Furthermore, although Rui et al. (
      • Rui H.
      • Artigas P.
      • Roux B.
      The selectivity of the Na+/K+-pump is controlled by binding site protonation and self-correcting occlusion.
      ) 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 sodium- and 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 Asp926 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 Asp926 protonation by changing local electrostatic potential and increasing cytoplasmic accessibility of a cavity in the protein that extends to Asp926. Access to the cytoplasmic anion-binding 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.

      Experimental procedures

      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 (
      • Jo S.
      • Kim T.
      • Iyer V.G.
      • Im W.
      CHARMM-GUI: a web-based graphical user interface for CHARMM.
      ) 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.
      Table 1Summary of the system specifics for molecular dynamics simulations
      Total number of atomsNumber of ionsNumber of lipid moleculesNumber of water molecules
      Na+-E1P233,660Na+, 185; Cl, 16431456,953
      K+-E2P261,581K+, 207; Cl, 18533665,179

      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 (
      • Essmann U.
      • Perera L.
      • Berkowitz M.
      • Darden T.
      • Lee H.
      • Pedersen L.
      A smooth particle mesh Ewald method.
      ) 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 (
      • Ryckaert J.
      • Ciccotti G.
      • Berendsen H.
      Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.
      ). 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 (
      • Phillips J.C.
      • Braun R.
      • Wang W.
      • Gumbart J.
      • Tajkhorshid E.
      • Villa E.
      • Chipot C.
      • Skeel R.D.
      • Kalé L.
      • Schulten K.
      Scalable molecular dynamics with NAMD.
      ). The reversible reference system propagator algorithms (rRESPA) multiple time step method was used (
      • Tuckerman M.
      • Berne B.J.
      • Martyna G.J.
      Reversible multiple time scale molecular dynamics.
      ) with a high frequency time step of 2.0 fs and a low frequency time step of 4.0 fs. Molecular configurations were saved every 20 ps. Harmonic restraints with a force constant of 5 kcal/mol/Å2 were applied on the phosphorus and nitrogen atoms of lipids and protein backbones and released gradually during the first 10 ns of simulations. Simulation times for equilibration runs are shown in supplemental 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, Asp926 was protonated in the E2P state to compare the results with the deprotonated Asp926 case (see Fig. 1 and the main text). In all cases, two distanceZ collective variables (
      • Fiorin G.
      • Klein M.L.
      • Hénin J.
      Using collective variables to drive molecular dynamics simulations.
      ) 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 − λ)Uprot + λUdeprot, 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 pKa 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 Asp804 and Glu779 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 reaction procedure with this alternative approach for the two FEP reactions in the potassium-bound E2P state. We chose Asp269 from the β-subunit of the Na+/K+-ATPase to protonate/deprotonate during FEP reactions to cancel charge-changing effects in simulations where the charge of the binding site changes. Asp269 from β-subunit (hereafter denoted as Asp269_B) is far from the binding site and in the extracellular side of the membrane; therefore, the effect of this residue on the binding site is most likely miniscule. The results of these FEP reactions are presented in supplemental Table S10. In the case of the Glu327-Asp269_B → Glu327-Asp808 transformation, the FEP result from this alternative approach agrees qualitatively with our reference reaction FEP approach. In the other case, the Glu327-Glu779-Asp269_B → Glu327-Glu779-Asp808 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 Glu327-Glu779 and Glu327-Glu779-Asp808 protonation states have similar free energies, and both have higher energy than Glu327-Asp808. 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 sodium-bound 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(AB) = −kBT ln〈exp(−(EBEA)/kBT)〉A (
      • Zwanzig R.W.
      High-temperature equation of state by a perturbation method. I. Nonpolar gases.
      ). 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 (
      • Grossfield A.
      • Ren P.
      • Ponder J.W.
      Ion solvation thermodynamics from simulation with a polarizable force field.
      )). 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.
      E1P3Na+ΔG(E1P)E1P3K+3Na+(water)ΔG(solution)3K+(water)REACTION 1


      Selectivities are reported as ΔΔG = ΔG(E1P) − ΔG(solution).

      Analysis of FEP simulations

      The multistate Bennett acceptance ratio method was used for free-energy estimation as implemented in the pymbar Python package (
      • Klimovich P.V.
      • Shirts M.R.
      • Mobley D.L.
      Guidelines for the analysis of free energy calculations.
      ,
      • Shirts M.R.
      • Chodera J.D.
      Statistically optimal analysis of samples from multiple equilibrium states.
      ). 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 Asp926, Glu327, Glu779, Asp804, and Asp808, 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.
      U=Ωdrρ(r)Vext(r)t
      (Eq. 1)


      where ρ(r) is the density of the atoms belonging to side chains of residues Asp926, Glu327, Glu779, Asp804, or Asp808; 〈Vext(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 (
      • Essmann U.
      • Perera L.
      • Berkowitz M.
      • Darden T.
      • Lee H.
      • Pedersen L.
      A smooth particle mesh Ewald method.
      ) 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, 〈Vext(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.
      ρ(r)=i1σi3(2π)32e(r-ri(t))22σ2t
      (Eq. 2)


      In this equation, r indicates the position in space, ri(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 (
      • Schmidtke P.
      • Bidon-Chanal A.
      • Luque F.J.
      • Barril X.
      MDpocket: open-source cavity detection and characterization on molecular dynamics trajectories.
      ,
      • Le Guilloux V.
      • Schmidtke P.
      • Tuffery P.
      Fpocket: an open source platform for ligand pocket detection.
      ), which utilizes Voronoi tessellation to define protein volumes. We used the tool MDpocket to analyze the structure of the protein along the trajectory and identified transient and stable pockets (
      • Schmidtke P.
      • Bidon-Chanal A.
      • Luque F.J.
      • Barril X.
      MDpocket: open-source cavity detection and characterization on molecular dynamics trajectories.
      ).

      Author contributions

      A. M. R. conducted most of the simulations, analyzed data, and wrote the initial draft. V. C. and L. D. designed and supervised the simulation protocols. A. M. R., V. A. V., and V. C. conceived the project and experimental design. J. R. B. helped in integrating simulation results with experimental observations. All authors participated in interpreting the data and in writing the manuscript.

      Supplementary Material

      References

        • Skou J.C.
        The identification of the sodium pump.
        Biosci. Rep. 1998; 18: 155-169
        • Skou J.C.
        The influence of some cations on an adenosine triphosphatase from peripheral nerves.
        Biochim. Biophys. Acta. 1957; 23: 394-401
        • Jorgensen P.L.
        • Hakansson K.O.
        • Karlish S.J.
        Structure and mechanism of Na,K-ATPase: functional sites and their interactions.
        Annu. Rev. Physiol. 2003; 65: 817-849
        • Schwinger R.H.
        • Bundgaard H.
        • Müller-Ehmsen J.
        • Kjeldsen K.
        The Na, K-ATPase in the failing human heart.
        Cardiovasc. Res. 2003; 57: 913-920
        • Poulsen H.
        • Khandelia H.
        • Morth J.P.
        • Bublitz M.
        • Mouritsen O.G.
        • Egebjerg J.
        • Nissen P.
        Neurological disease mutations compromise a C-terminal ion pathway in the Na+/K+-ATPase.
        Nature. 2010; 467: 99-102
        • Moseley A.E.
        • Williams M.T.
        • Schaefer T.L.
        • Bohanan C.S.
        • Neumann J.C.
        • Behbehani M.M.
        • Vorhees C.V.
        • Lingrel J.B.
        Deficiency in Na,K-ATPase α isoform genes alters spatial learning, motor activity, and anxiety in mice.
        J. Neurosci. 2007; 27: 616-626
        • Morth J.P.
        • Poulsen H.
        • Toustrup-Jensen M.S.
        • Schack V.R.
        • Egebjerg J.
        • Andersen J.P.
        • Vilsen B.
        • Nissen P.
        The structure of the Na+, K+-ATPase and mapping of isoform differences and disease-related mutations.
        Philos. Trans. R. Soc. Lond. B Biol. Sci. 2009; 364: 217-227
        • Yatime L.
        • Buch-Pedersen M.J.
        • Musgaard M.
        • Morth J.P.
        • Lund Winther A.M.
        • Pedersen B.P.
        • Olesen C.
        • Andersen J.P.
        • Vilsen B.
        • Schiøtt B.
        • Palmgren M.G.
        • Møller J.V.
        • Nissen P.
        • Fedosova N.
        P-type ATPases as drug targets: tools for medicine and science.
        Biochim. Biophys. Acta. 2009; 1787: 207-220
        • Post R.L.
        • Sen A.K.
        • Rosenthal A.S.
        A phosphorylated intermediate in adenosine triphosphate-dependent sodium and potassium transport across kidney membranes.
        J. Biol. Chem. 1965; 240: 1437-1445
        • Albers R.
        Biochemical aspects of active transport.
        Annu. Rev. Biochem. 1967; 36: 727-756
        • Post R.L.
        • Kume S.
        • Tobin T.
        • Orcutt B.
        • Sen A.K.
        Flexibility of an active center in sodium-plus-potassium adenosine triphosphatase.
        J. Gen. Physiol. 1969; 54: 306-326
        • Kanai R.
        • Ogawa H.
        • Vilsen B.
        • Cornelius F.
        • Toyoshima C.
        Crystal structure of a Na+-bound Na+,K+-ATPase preceding the E1P state.
        Nature. 2013; 502: 201-206
        • Nyblom M.
        • Poulsen H.
        • Gourdon P.
        • Reinhard L.
        • Andersson M.
        • Lindahl E.
        • Fedosova N.
        • Nissen P.
        Crystal structure of Na+, K+-ATPase in the Na+-bound state.
        Science. 2013; 342: 123-127
        • Laursen M.
        • Yatime L.
        • Nissen P.
        • Fedosova N.U.
        Crystal structure of the high-affinity Na+,K+-ATPase–ouabain complex with Mg2+ bound in the cation binding site.
        Proc. Natl. Acad. Sci. U.S.A. 2013; 110: 10958-10963
        • Morth J.P.
        • Pedersen B.P.
        • Toustrup-Jensen M.S.
        • Sørensen T.L.
        • Petersen J.
        • Andersen J.P.
        • Vilsen B.
        • Nissen P.
        Crystal structure of the sodium-potassium pump.
        Nature. 2007; 450: 1043-1049
        • Ogawa H.
        • Shinoda T.
        • Cornelius F.
        • Toyoshima C.
        Crystal structure of the sodium-potassium pump (Na+,K+-ATPase) with bound potassium and ouabain.
        Proc. Natl. Acad. Sci. U.S.A. 2009; 106: 13742-13747
        • Shinoda T.
        • Ogawa H.
        • Cornelius F.
        • Toyoshima C.
        Crystal structure of the sodium-potassium pump at 2.4 Å resolution.
        Nature. 2009; 459: 446-450
        • Møller J.
        • Nissen P.
        • Sørensen T.
        • le Maire M.
        Transport mechanism of the sarcoplasmic reticulum Ca2+-ATPase pump.
        Curr. Opin. Struct. Biol. 2005; 15: 387-393
        • Yu H.
        • Ratheal I.M.
        • Artigas P.
        • Roux B.
        Protonation of key acidic residues is critical for the K+-selectivity of the Na/K pump.
        Nat. Struct. Mol. Biol. 2011; 18: 1159-1163
        • Milanick M.A.
        • Arnett K.L.
        Extracellular protons regulate the extracellular cation selectivity of the sodium pump.
        J. Gen. Physiol. 2002; 120: 497-508
        • Li H.
        • Ngo V.
        • Da Silva M.C.
        • Salahub D.R.
        • Callahan K.
        • Roux B.
        • Noskov S.Y.
        Representation of ion–protein interactions using the Drude polarizable force-field.
        J. Phys. Chem. B. 2015; 119: 9401-9416
        • Li H.
        • Robertson A.D.
        • Jensen J.H.
        Very fast empirical prediction and rationalization of protein pKa values.
        Proteins. 2005; 61: 704-721
        • Fiorin G.
        • Klein M.L.
        • Hénin J.
        Using collective variables to drive molecular dynamics simulations.
        Mol. Phys. 2013; 111: 3345-3362
        • Pedersen P.A.
        • Nielsen J.M.
        • Rasmussen J.H.
        • Jorgensen P.L.
        Contribution to Tl+, K+, and Na+ binding of Asn776, Ser775, Thr774, Thr772, and Tyr771 in cytoplasmic part of fifth transmembrane segment in R-subunit of renal Na,K-ATPase.
        Biochemistry. 1998; 37: 17818-17827
        • Klodos I.
        • Post R.L.
        • Forbush 3rd, B.
        Kinetic heterogeneity of phosphoenzyme of Na,K-ATPase modeled by unmixed lipid phases.
        J. Biol. Chem. 1994; 269: 1734-1743
        • Mahmmoud Y.A.
        • Kopec W.
        • Khandelia H.
        K+ congeners that do not compromise Na+ activation of the Na+, K+-ATPase: hydration of the ion binding cavity likely controls ion selectivity.
        J. Biol. Chem. 2015; 290: 3720-3731
        • Jorgensen P.L.
        • Rasmussen J.H.
        • Nielsen J.M.
        • Pedersen P.A.
        Transport-linked conformational changes in Na,K-ATPase. Structure-function relationships of ligand binding and E1-E2 conformational transitions.
        Ann. N.Y. Acad. Sci. 1997; 834: 161-174
        • Lim H.H.
        • Shane T.
        • Miller C.
        Intracellular proton access in a Cl/H+ antiporter.
        PLoS Biol. 2012; 10: e1001441
        • Ogawa H.
        • Qiu Y.
        • Philo J.S.
        • Arakawa T.
        • Ogata C.M.
        • Misono K.S.
        Reversibly bound chloride in the atrial natriuretic peptide receptor hormone-binding domain: possible allosteric regulation and a conserved structural motif for the chloride-binding site.
        Protein Sci. 2010; 19: 544-557
        • Li C.
        • Geering K.
        • Horisberger J.
        The third sodium binding site of Na,K-ATPase is functionally linked to acidic pH-activated inward current.
        J. Membr. Biol. 2006; 213: 1-9
        • Yaragatupalli S.
        • Olivera J.F.
        • Gatto C.
        • Artigas P.
        Altered Na+ transport after an intracellular α-subunit deletion reveals strict external sequential release of Na+ from the Na/K pump.
        Proc. Natl. Acad. Sci. U.S.A. 2009; 106: 15507-15512
        • Riant F.
        • De Fusco M.
        • Aridon P.
        • Ducros A.
        • Ploton C.
        • Marchelli F.
        • Maciazek J.
        • Bousser M.G.
        • Casari G.
        • Tournier-Lasserve E.
        ATP1A2 mutations in 11 families with familial hemiplegic migraine.
        Hum. Mutat. 2005; 26: 281
        • Schmidtke P.
        • Bidon-Chanal A.
        • Luque F.J.
        • Barril X.
        MDpocket: open-source cavity detection and characterization on molecular dynamics trajectories.
        Bioinformatics. 2011; 27: 3276-3285
        • Tavraz N.N.
        • Friedrich T.
        • Dürr K.L.
        • Koenderink J.B.
        • Bamberg E.
        • Freilinger T.
        • Dichgans M.
        Diverse functional consequences of mutations in the Na+/K+-ATPase α2-subunit causing familial hemiplegic migraine type 2.
        J. Biol. Chem. 2008; 283: 31097-31106
        • Weigand K.M.
        • Swarts H.G.
        • Russel F.G.
        • Koenderink J.B.
        Biochemical characterization of sporadic/familial hemiplegic migraine mutations.
        Biochim. Biophys. Acta. 2014; 1838: 1693-1700
        • Toustrup-Jensen M.S.
        • Holm R.
        • Einholm A.P.
        • Schack V.R.
        • Morth J.P.
        • Nissen P.
        • Andersen J.P.
        • Vilsen B.
        The C terminus of Na+,K+-ATPase controls Na+ affinity on both sides of the membrane through Arg935.
        J. Biol. Chem. 2009; 284: 18715-18725
        • Suzuki K.
        • Post R.
        Equilibrium of phosphointermediates of sodium and potassium ion transport adenosine triphosphatase. Action of sodium ion and Hofmeister effect.
        J. Gen. Physiol. 1997; 109: 537-554
        • Ganea C.
        • Babes A.
        • Lüpfert C.
        • Grell E.
        • Fendler K.
        • Clarke R.J.
        Hofmeister effects of anions on the kinetics of partial reactions of the Na+,K+-ATPase.
        Biophys. J. 1999; 77: 267-281
        • Tadini-Buoninsegni F.
        • Moncelli M.R.
        • Peruzzi N.
        • Ninham B.W.
        • Dei L.
        • Nostro P.L.
        Hofmeister effect of anions on calcium translocation by sarcoplasmic reticulum Ca2+-ATPase.
        Sci. Rep. 2015; 5: 14282
        • Clarke R.J.
        Dipole-potential-mediated effects on ion pump kinetics.
        Biophys. J. 2015; 109: 1513-1520
        • Rui H.
        • Artigas P.
        • Roux B.
        The selectivity of the Na+/K+-pump is controlled by binding site protonation and self-correcting occlusion.
        Elife. 2016; 5: e16616
        • Jo S.
        • Kim T.
        • Iyer V.G.
        • Im W.
        CHARMM-GUI: a web-based graphical user interface for CHARMM.
        J. Comput. Chem. 2008; 29: 1859-1865
        • Essmann U.
        • Perera L.
        • Berkowitz M.
        • Darden T.
        • Lee H.
        • Pedersen L.
        A smooth particle mesh Ewald method.
        J. Chem. Phys. 1995; 103: 8577
        • Ryckaert J.
        • Ciccotti G.
        • Berendsen H.
        Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.
        J. Comput. Phys. 1977; 23: 327-341
        • Phillips J.C.
        • Braun R.
        • Wang W.
        • Gumbart J.
        • Tajkhorshid E.
        • Villa E.
        • Chipot C.
        • Skeel R.D.
        • Kalé L.
        • Schulten K.
        Scalable molecular dynamics with NAMD.
        J. Comput. Chem. 2005; 26: 1781-1802
        • Tuckerman M.
        • Berne B.J.
        • Martyna G.J.
        Reversible multiple time scale molecular dynamics.
        J. Chem. Phys. 1992; 97: 1990-2001
        • Zwanzig R.W.
        High-temperature equation of state by a perturbation method. I. Nonpolar gases.
        J. Chem. Phys. 1954; 22: 1420
        • Grossfield A.
        • Ren P.
        • Ponder J.W.
        Ion solvation thermodynamics from simulation with a polarizable force field.
        J. Am. Chem. Soc. 2003; 125: 15671-15682
        • Klimovich P.V.
        • Shirts M.R.
        • Mobley D.L.
        Guidelines for the analysis of free energy calculations.
        J. Comput. Aided Mol. Des. 2015; 29: 397-411
        • Shirts M.R.
        • Chodera J.D.
        Statistically optimal analysis of samples from multiple equilibrium states.
        J. Chem. Phys. 2008; 129: 124105
        • Le Guilloux V.
        • Schmidtke P.
        • Tuffery P.
        Fpocket: an open source platform for ligand pocket detection.
        BMC Bioinformatics. 2009; 10: 168