The Na+/K+-pump maintains the physiological K+ and Na+ electrochemical gradients across the cell membrane. It operates via an 'alternating-access' mechanism, making iterative transitions between inward-facing (E1) and outward-facing (E2) conformations. Although the general features of the transport cycle are known, the detailed physicochemical factors governing the binding site selectivity remain mysterious. Free energy molecular dynamics simulations show that the ion binding sites switch their binding specificity in E1 and E2. This is accompanied by small structural arrangements and changes in protonation states of the coordinating residues. Additional computations on structural models of the intermediate states along the conformational transition pathway reveal that the free energy barrier toward the occlusion step is considerably increased when the wrong type of ion is loaded into the binding pocket, prohibiting the pump cycle from proceeding forward. This self-correcting mechanism strengthens the overall transport selectivity and protects the stoichiometry of the pump cycle.https://doi.org/10.7554/eLife.16616.001
A protein called the sodium-potassium pump resides in the membrane that surrounds living cells. The role of this protein is to 'pump' sodium and potassium ions across the membrane to help restore their concentration inside and outside of the cell. About 25% of the body's energy is used to keep the pump going, rising to nearly 70% in the brain. Problems that affect the pump have been linked to several disorders, including heart, kidney and metabolic diseases, as well as severe neurological conditions.
The sodium-potassium pump must be able to effectively pick out the correct ions to transport from a mixture of many different types of ions. However, it was not clear how the pump succeeds in doing this efficiently. Rui et al. have now used a computational method called molecular dynamics simulations to model how the sodium-potassium pump transports the desired ions across the cell membrane.
The pump works via a so-called 'alternating-access' mechanism, repeatedly transitioning between inward-facing and outward-facing conformations. In each cycle, it binds three sodium ions from the cell’s interior and exports them to the outside. Then, the pump binds to two potassium ions from outside the cell and imports them inside. Although the bound sodium and potassium ions interact with similar binding sites in the pump, the pump sometimes preferentially binds sodium, and sometimes potassium. The study performed by Rui et al. shows that this preference is driven by how protons (hydrogen ions) bind to the amino acids that make up the binding site.
The simulations also suggest that the pump uses a ‘self-correcting’ mechanism to prevent the pump from transporting the wrong types of ions. When incorrect ions are present at the binding sites, the pump cycle pauses temporarily until the ions detach from the pump. Only when the correct ions are bound will the pump cycle continue again.
In the future, Rui et al. hope to use long time-scale molecular dynamics simulations to show the conformational transition in action. In addition, the 'self-correcting' mechanism will be directly tested by letting the wrong and correct ions compete for the binding sites to see whether the pump will transport only the correct ions.https://doi.org/10.7554/eLife.16616.002
The Na+/K+-pump is a primary active membrane transporter present in nearly all animal cells. It belongs to the P-type ATPase family, which utilizes the energy released from ATP hydrolysis to move ions against their concentration gradients across a membrane barrier. The ion species transported by the pump under physiological conditions are Na+ and K+. Like many other membrane transporters, the Na+/K+-pump works according to an 'alternating-access' ion transport mechanism, with the bound ions accessible from only one side of the membrane at a time. The consensus scheme of the pump cycle is known as the 'Albers-Post' cycle (Albers, 1967; Post et al., 1969). It involves two major conformations, E1 and E2, with inward- and outward-facing ion binding sites, respectively. In each cycle, the E1 conformation binds three Na+ from the cytosol and exports them using the energy from ATP hydrolysis. After external release of Na+, binding of extracellular K+ promotes the structural transition to the occluded E2 state, which finally imports two K+ as binding of ATP returns the pump to the E1 conformation. E1/E2 conformational change and the accompanied electrogenic active transport are facilitated by the phosphorylation and dephosphorylation of an aspartate residue in the cytoplasmic domain, which is a hallmark of the P-type ATPase family (Axelsen et al., 1998).
The presence of the Na+/K+-pump is essential for excitability and secondary active transport. More than 40% of the energy produced in mammals is consumed by the Na+/K+-pump (Milligan et al., 1985). Although it is a machine designed for such a precise and important function, it has been shown that many cations, including congeners of K+ and some organic cations, can bind to the same sites used by the pump to bind and transport K+ ions (Mahmmoud et al., 2015; Ratheal et al., 2010). Competitive binding between Na+ and other cations at the cytoplasmic side of the membrane has also been observed (Schneeberger et al., 2001). An unsolved puzzle, therefore, is how the Na+/K+-pump is able to recognize and accept two K+ from the extracellular matrix, where Na+ concentration is much higher, and how it selectively binds Na+ from the cytoplasm to keep the pump cycle running forward.
Structural studies have provided tremendous insights into the function of the Na+/K+-pump, which consists of two obligatory subunits, α (catalytic) and β (auxiliary), and sometimes a tissue specific regulatory subunit from the FXYD protein family (Geering, 2006; Mercer et al., 1993). The transmembrane region of the α-subunit contains the ion binding sites within its 10 helices (called M1-M10). Recent crystal structures of the pump in its E1 and E2 states reveal the locations of the three Na+ binding sites in E1 and the two K+ binding sites in E2 (Kanai et al., 2013; Laursen et al., 2013; Morth et al., 2007; Shinoda et al., 2009). From structural alignments based on the heavy atom positions in the binding sites, it becomes clear that the binding pocket harboring sites I and II in E1 overlaps with those in E2 (Figure 1). Site III is only formed in E1. It is located between the transmembrane helices M5, M6, and M8 and is thought to exclusively bind Na+ but appears to catalyze H+ transport (Poulsen et al., 2010; Ratheal et al., 2010) in a manner that presents a complex dependence on the concentrations of Na+, K+, and H+ (Mitchell et al., 2014). Previous studies have indicated that protonation at the ion binding pocket may play a role in the selectivity of these sites for K+ when the pump is in its E2 state (Yu et al., 2011). Biochemical assays on the mutant D926N, which is often used as a surrogate for protonated D926, also show that it induces distinct effects on Na+ and K+ binding (Einholm et al., 2010; Jewell-Motz et al., 1993), suggesting a potential change in its protonation state upon the E1/E2 transition. Taken together, the evidence, although indirect, suggests the possibility of an E1 specific protonation state that favors Na+ binding to the pump. The nature of this protonation state is, however, unclear.
In the present study, we started with the recently published crystal structure of the Na+/K+-pump in a partially occluded Na3·E1(ADP·Pi) state and used molecular dynamics (MD) simulations to show that there is a correlation between the binding pocket protonation and the Na+ selectivity in E1. The binding sites in Na3·E1(ADP·Pi) were tested one by one by 'alchemically transforming' the bound Na+ into a K+ while keeping the other two sites occupied by Na+. A similar approach was previously used to study the selectivity in other conformational states of the pump (Yu et al., 2011) (see also Materials and methods section). The results show that the Na+ selectivity at all three sites is realized only when a specific set of binding pocket residues are protonated. This set of residues is, by comparison, different from that in the K+ selective E2 state. The implication is that the protonation state and the selectivity of the pump are tightly coupled; when the pump undergoes a transition between E1 and E2, a protonation state switch occurs. The present findings also show that the effective selectivity of the pump is reinforced by a self-correcting mechanism, which prevents the occlusion step on either side of the membrane if the wrong type of ions were loaded into the binding sites. This mechanism ensures that counterproductive transport cycles do not occur.
The crystal structures of the pump in its E1 (PDBID 3WGV [Kanai et al., 2013]) and E2 (PDBID 2ZXE [Shinoda et al., 2009]) conformations show that there is a large structural overlap between the Na+ and K+ binding pockets. The sites I and II are in the main binding chamber formed by helices M4, M5, and M6, while site III, located in between M5, M6, and M8 (Figure 1), is Na+ exclusive and appears in E1 only. The coordination of Na+ and K+ at sites I and II are similar. Many residues are found to coordinate both Na+ and K+ at these two sites in E1 and E2 (Table 1). Although the structural difference between the ion binding pockets in E1 and E2 may not be particularly large, the local physicochemical environment can display considrable variations. In particular, the latter could affect the pKa and the protonation states of the six protonatable residues in the binding pocket (Figure 1). Their pKa values calculated using PROPKA3.1 (Olsson et al., 2011) are listed in Table 2.
It is important to note that the pKa values calculated with empirical methods like PROPKA are sensitive to local structural perturbations. Even when the structures assume the identical conformational state and are taken from the same asymmetric unit from a crystal, the pKa values of the same residue can differ by more than one pH unit. For example, the pKa of E779 is 9.9 in chain A of the crystal structure 3WGV but the value is 8.4 in chain C, yet, the structural difference between the two chains is minimal (Table 2). Because of this, only the structures with the highest resolution for the E1 and E2 states of the pump, 3WGV and 2ZXE, were used to guide the protonation state assignments in order to avoid false assignments of protonation state originated from structural uncertainty in lower resolution crystal structures. Interestingly, the protonation states of E327 and D804 assigned based on the 3WGV and 4HQJ structures are reversed. According to the PROPKA results, E327 appears to be deprotonated and D804 protonated in 3WGV, and vice versa in 4HQJ. While this inconsistency could be due to the lower resolution in the crystal structure of 4HQJ, it is worth noting that these two residues are in close proximity from one another suggesting that a proton could be passed back and forth between them in the E1 state.
Previous calculations have indicated that a specific set of residues has to be protonated for the E2 binding sites to be K+ selective (Yu et al., 2011). Although the protonation states of D926 and E954 were not considered in the previous study because they lie outside of the two K+ binding sites, they are likely to be protonated in E2 according to the predicted pKa (Table 2). The binding pocket protonation derived from the crystal structure of the pump in its partially occluded Na3·E1(ADP·Pi) state differs from that in E2. According to their pKa values, D804, D808, and D926 should be deprotonated and the glutamates, E327, E786, and E954, should be protonated. Arguments based on pKa values predicted by PROPKA, however, have to be taken with caution, because the empirical method is highly sensitive to the local structural variation. For example, a deprotonated acidic side chain could have a PROPKA predicted pKa value at slightly higher than 7 because of the structural snapshot used to make the prediction. To seek a more robust assessment of these factors, molecular dynamics (MD) simulations were conducted to examine the structural stability of the binding pocket for different protonation configurations. One goal is to determine the protonation states of D804 and D926, both of which have a predicted pKa value within 1.5 pH unit to 7. The protonation state of D808 was also investigated since it affects the K+ selectivity of the binding sites in the E2(K2) state. A total of eight MD simulations were carried out (Table 3), starting from all possible protonation state configurations accessible by these three residues. The protocol for setting up these systems is given in details in the Materials and methods section.
Figure 2 shows the binding pocket conformation in the E1 systems (Table 3) at the end of the all-atom MD simulations. Out of the three aspartates when two or more are protonated the binding pocket becomes unstable, showing a large deviation from the crystal structure either with one of the bound Na+ expelled from its binding site (systems E1_S5-7), or with a Cl- ion attracted into the binding pocket (system E1_S4). These scenarios are not likely to happen in the normal function cycle of the pump. On the other hand, when the number of protonated aspartates is less than two the binding pocket remains structurally close to that in the crystal structure (systems E1_S0-3). The heavy atom root mean squared deviations (RMSD) between the crystal structure and the structure snapshots from the last 50 ns of the simulations are less than 1.6 Å.
Analysis of the structural stability of the binding sites based on the MD simulations indicates that all of the four protonation states producing stable ion binding pockets may coexist when the pump is in the Na3·E1(ADP·Pi) state. The relative population of these states in the state Na3·E1(ADP·Pi), however, would differ from one another. The overall selectivity of the binding sites has contributions from all the states, and the predominant protonation state configuration should produce binding sites that are Na+ selective. The determination of this protonation state configuration starts from the protein-membrane systems generated by the restraint-free MD simulations. A reduced structural model of the binding pocket is derived from each of these systems (Table 3) and the binding free energy differences between Na+ and K+ () at the binding sites are calculated (See Materials and methods). The value of reflects the affinity difference between K+ and Na+ binding as . The results are plotted as the logarithm of the affinity ratio, , in Figure 3. The values of are shown in Table 4.
As shown in Figure 3, the protonation state of the aspartates greatly affects both the E1 and E2 binding sites selectivity. The protonation state at E2(K2) binding pocket proposed in the earlier study (Yu et al., 2011) is confirmed by the present calculations, in which the reduced system includes a larger region of the protein and more residues at the outskirt of the main binding pocket, including D926 and E954. It is evident that a different set of residues is protonated in the Na+ selective E1 as compared to in the K+ selective E2 (Figure 3). While the glutamates (i.e., E327, E779, and E954) remain protonated in both the Na3·E1(ADP·Pi) and E2(K2) states, the binding pocket aspartates take on opposite protonation states. In Na3·E1(ADP·Pi) D804 is protonated and both D808 and D926 are deprotonated, and their protonation states are reversed in E2(K2). Free energy perturbation (FEP) calculations on the protonation/deprotonation of D804 estimated a pKa value of 9.2 (See Materials and methods), further supporting its protonated form under physiological conditions in E1.
The impact of charge-neutralizing mutations of the key aspartate residues on the binding site selectivity was examined. D to N mutations have previously been used in experimental studies as a strategy to ascertain the possible charged state of these residues. The results of the free energy calculations, summarized in Figure 4, indicate that all the three sites remain Na+ selective in the partly occluded Na3·E1(ADP·Pi) in both D804N and D808N, whilst D926N causes the sites to lose most of the Na+ selectivity. Using Na+ dependent phosphorylation and ATP binding assays, it was shown that the cytoplasmic binding affinity for both Na+ and K+ is reduced in D804N and D808N (Jorgensen et al., 2001; Pedersen et al., 1997). The D804N mutation affects the cytoplasmic K+ and Na+ binding differently. The mutation’s impact on cytoplasmic Na+ binding is not as prominent as on K+ binding and D804N would appear more Na+ selective than the wildtype. This is reflected in the calculation results. The KD,K/KD,Na at sites I and III are reduced by a factor of ~50 from 476.6 and 2523.3 in the wildtype to 10.3 and 54.6 in the D804N mutant, but the sites are still Na+ selective. The selectivity of site II is increased to a much larger extent, from 17.0 in the wildtype to 13359.7 in the mutant. This is an ~1000 fold increase in selectivity and it implies that D804N is likely more Na+ selective than the wildtype pump. The D808N mutant also becomes more Na+ selective compared to the wildtype pump. The most dramatic increase in Na+ selectivity happens at site II as in the case of D804N. Interestingly, site I, which prefers to bind K+ in the wildtype with the protonated D808 (Figure 3), is now Na+ selective in D808N. It could be that the spatial packing is the major contributing factor to the site I ion selectivity. In the calculations, both sites II and III lose their Na+ selectivity in the D926N mutant. The shifting of selectivity towards K+ at site III in this mutant is a bit surprising, and worth commenting on. There are two possible explanations. First, substituting a negatively charged carbonyl oxygen with a bulkier but neutral –NH2 at the D926 sidechain could prevent entrance of an ion to this site. Thus, the simulated conformation with an ion at this site could be energetically inaccessible. In other words, this is a metastable state with the absolute free energy of Na+ and K+ binding to this conformation equally prohibitive. An alternative explanation is that the D926N mutation alters the available conformational space accessible by helix M5 and this allows K+ to go into site III as suggested in reference (Kanai et al., 2013). The binding of this K+, however, prevents the further binding of Na+, and results in the compromised Na+ selectivity in this mutant. Even though the calculations show D926N with decreased Na+ selectivity, they do not explain why experimentally the D926N mutant has compromised Na+ binding ability in the absence of K+ (Einholm et al., 2010). The loss in selectivity seen in this D926N mutant, however, could account for the strong inhibition by high K+ on Na+/K+-ATPase activity (cf. Figure 2A in reference [Einholm et al., 2010]).
All the D to N mutations tested for the occluded E2(K2) state appear to have little impact on the K+ selectivity, in apparent contradiction with the experimental observations showing that both D804N and D808N display decreased selectivity for external K+ (Kuntzweiler et al., 1996; Pedersen et al., 1997). In a previous computational study, the D808N mutation was also shown to minimally affect the K+ selectivity in E2(K2). Discrepancies between the calculations and the selectivity inferred from biochemical experiments have been noted previously, though the reason was not clear. A plausible explanation could be that the experimentally measured selectivity is an outcome from the relevant states including both P-E2 and E2(K2). However, any contribution from the P-E2 state to the observed selectivity was not taken into account by the calculations because of the missing crystal structure of the outward-facing state. To test this idea, we generated a model structure of the P-E2 state based on the homolog structures from the Ca2+ SERCA pump (see below). This model was then used to calculate the in the wildtype and the mutant pumps. The results show that both sites I and II have increased preference for external Na+ in the P-E2 state of D808N (Figure 4), which explains previous discrepancies between the calculations and the experiments.
Using the crystal structures of the Ca2+ SERCA pump as templates and a coarse grained transition pathway calculation method, ANMpathway (Das et al., 2014), we generated structural models of the Na+– and K+–loaded outward-facing P-E2 state (see Materials and methods). The models appear to be structurally similar to the recently published P-E2 like structures of the Na+/K+-pump (Gregersen et al., 2016; Laursen et al., 2013) and the vanadate-Inhibited, P-E2 mimic of the Ca2+ SERCA pump (Clausen et al., 2016). The MD simulations of the models are also able to predict with considerable accuracy the experimentally measured gating charge upon ion binding (Castillo et al., 2015). The structural transition pathways leading to these models provide a view of the intermediate structures along the pump cycle. Using these models, we calculated the for the binding site of the intermediate state structures. The entire atomistic protein-membrane systems were used in the calculations. The results are presented in Figure 5 and Table 5. The calculations are valid as the values mirror those computed from the reduced systems with the same protonation states (Figure 3).
The results reveal a fascinating feature of the binding site selectivity along the pump cycle. Initially, the two K+ binding sites in the outward-facing model P-E2·K2 (state ① in Figure 5) appear to be non-selective or only weakly K+ selective. However, the selectivity for K+ over Na+ increases as the pump occludes to form the K+–bound occluded state E2(K2) (state ③). After the dephosphorylation of P-E2, the binding pocket opens up to the cytoplasmic side. Accompanying this structural transition is the selectivity reversal at site I, switching the site from K+ to Na+ selective (③ to ④ in Figure 5). The protonation state change further shifts the ion selectivity at the sites in E1 towards Na+. Changing the protonation state in K2E1 (④ to ⑤ in Figure 5) to the protonation state dominant in Na3·E1(ADP·Pi), i.e., state ⑥ in Figure 5 with protonated D804 and deprotonated D808 and D926, further reduces the K+ selectivity at both sites I and II. When the pump is in this state, site I is Na+ selective and site II is only weakly K+ selective. Together, the E2 to E1 structural transition and the protonation state change promote the release of K+ to the cytoplasmic side. Transient changes in the binding pocket protonation state upon occlusion/deocclusion are possible and could lead to variations in the magnitude of for the intermediate states. Nonetheless, the general trend in free energy changes upon occlusion/deocclusion should remain because the of the occluded state E2(K2) has a much larger magnitude than that of the open-access outward- and inward-facing states.
The free energy calculations corroborate the notion that the two K+ binding from the extracellular side is sequential and possibly cooperative. The sequential binding of extracellular K+ was initially demonstrated by Forbush (Forbush, 1987). Recently using crystallography and isotopic measurements Ogawa et al. presented strong evidence that the first K+ binds at site I and the second K+ at site II (Ogawa et al., 2015). Whether there is cooperativity upon extracellular K+ binding, however, is not clear from the crystal structures. Two electrode voltage clamp experiments have shown that the two extracellular K+ binding events are relatively independent in the absence of extracellular Na+, but there is positive cooperativity of K+ binding when extracellular Na+ ions are present (Jaisser et al., 1994). The current calculations of the two ion bound P-E2 (Table 5) provide an interpretation of this phenomenon. With a K+ ion at site I, site II is relatively nondiscriminatory ( = 0.7 kcal/mol, p/p to p/s in state ① P-E2.K2, Table 5), but it becomes much more Na+ selective ( = 2.4 kcal/mol, s/s to s/p in state ⑨ P-E2·Na3, Table 5) when there is a Na+ ion bound at site I. In the presence of extracellular Na+, K+ ions have to first compete with Na+ to bind at site I in order for the subsequent K+ to bind, therefore resulting in the observed binding cooperativity in the presence of extracellular Na+.
A similar phenomenon is seen in the Na+ branch of the pump cycle. Initially, site I in the inward-facing Na3·E1(ADP·Pi) structure (state ⑦) is weakly Na+ selective when the other two sites are filled with Na+. However, the intermediate state ⑧ during the occlusion process shows that its site I is highly discriminating against K+ with a of 2.9 kcal/mol (Figure 5 and Table 5). Therefore, if a K+ is bound in this site, the free energy barrier toward occlusion would be much higher than when a Na+ ion was bound and the subsequent occlusion is not likely. In the case when a K+ ion replaces a Na+ at site II or III, although the energy barrier of occlusion is not as prominent, such a binding pocket ion configuration is energetically less favorable than the native configuration and its appearance is unlikely in the first place.
Unlike the K+ branch, the calculations along the Na+ branch offers limited insights on the sequence of Na+ binding from the cytoplasmic side. The results are indicative of the selectivity at the sites, but not the absolute binding affinity. Based on the calculations alone, it is not possible to determine which site is the first to bind cytoplasmic Na+. It could be site III, as suggested by Kanai et al., and this prepares the other two sites for the following Na+ binding (Kanai et al., 2013). Or, alternatively, the first two Na+ bind at sites I and II in the main binding pocket and the last Na+ enters, takes over site II, and pushes the two previously bound ions further, so that the ions in sites II and I now move to sites I and III, in a process reminiscent of the “knock-on” mechanism occurring in potassium channels (Hodgkin et al., 1955). A direct assessment of these proposed binding sequences will require further experiments.
The results from the MD simulations support the notion that the protonation state of the binding pocket and its selectivity are closely related. Because the binding pocket in the Na+/K+-pump displays considerable flexibility, it is worth pausing to reflect on the possible mechanism that underlies this relation. While the selectivity of a very rigid binding site is first and foremost predetermined by its atomic geometry, the selectivity of a flexible binding site is strongly affected by local ion-ligand and ligand-ligand interactions. In such flexible systems, selectivity is controlled by the both number and the physicochemical properties of ion-coordinating ligands (Yu et al., 2010). For example, high-field ligands, such as deprotonated acidic side chains tend to favor Na+ binding and protonation revert those to low-field ligands, which tend to favor K+. Hence, in the Na+/K+-pump protonation is exploited to modulate selectivity by altering the electrostatic properties of several of residues in the binding pocket. This is also consistent with the results from our previous study on the ion selectivity in E2(K2), which concluded that changes in the electrostatic properties of the protonatable residues was the likely mechanism responsible for the K+ selectivity in the E2 state of the pump (Yu et al., 2011). Even though the local structural rearrangements at the binding sites are small upon the change in protonation state, it is enough to change the physical properties of the coordinating ligand, thus giving rise to discernable differences in the ion selectivity. The crystal structures of the Na+/K+-pump in its Na3.E1·(ADP·Pi) and E2(K2) states show similarity in their binding pocket configurations (Figure 1), including the coordination patterns of the bound ions at sites I and II (Table 1). The heavy atom RMSD between the binding pocket residues is 2.5 Å and the structural difference remains after hundreds of ns of simulations. Empirical pKa and FEP calculations based on the MD simulation equilibrated structures indicate that the binding pocket glutamates (i.e., E327, E779, and E954) are likely protonated in both Na3.E1·(ADP·Pi) and E2(K2), although previous mutagenesis experiments showed that charge-neutralizing mutation E327Q affects pump function, possibly by altering ion-pump interactions and the kinetics of the occlusion/deocclusion reactions along the pump cycle (Jorgensen et al., 2001; Kuntzweiler et al., 1995; Li et al., 2006; Nielsen et al., 1998). The calculations also suggest that the protonation states of D804, D808, and D926 are different in Na3·E1·(ADP·Pi) and E2(K2). Among the three aspartates only D804 is protonated in order for all three sites to stay Na+ selective in Na3·E1·(ADP·Pi). The E1 binding pocket devoid of ions carries a net charge of -2, consistent with that deduced from previous fluorescence studies (Schneeberger et al., 1999). Mutagenesis experiments are also consistent with the unlikely protonation of D808 and D926 when the pump is in E1 trying to bind cytoplasmic Na+ (Jewell-Motz et al., 1993; Pedersen et al., 1997). The protonation states of these three residues are reversed in E2(K2) with D804 deprotonated and the other two protonated, a result that is supported by a previous computational study (Yu et al., 2011). The different protonation states for E1 and E2 also agree well with the 'four-site' model proposed by Skou and Esmann more than 30 years ago with the K+-bound E2 state carrying two sidechain protons (H2EK2) and the Na+-bound E1 state carrying only one (HENa3) (Skou et al., 1980).
It is worth pointing out that the second proton (i.e., the proton on D926) in the K+ bound E2 state must come from and return to the same side of the membrane during the pump cycle so that the net charge moved per cycle by the pump is one. It seems unlikely that this proton could come from the extracellular side because altering extracellular pH would then protonate or deprotonate D926, causing major changes in Na+ binding affinity and the maximum pumping turnover rate. This is not observed over an external pH range 9.6 to 5.6 (Mitchell et al., 2014; Vasilyev et al., 2004; Yu et al., 2011). Therefore, it is more likely that the D926 proton comes from the cytoplasmic side. This is supported by MD simulations of the P-E2·Na2 and Na2·E1·(ADP·Pi) revealing the existence of aqueous pathways connecting the cytoplasm and D926 (Figure 6). One of the water pathways is located between the helices M5, M7, and M8 (Figure 6B), similar to the C-terminal proton pathway previously proposed (Poulsen et al., 2010). A proton could enter through this pathway to protonate D926 and then leave through the same pathway during the E2 to E1 transition, or through an alternative path passing the main binding pocket along with the dissociating ions as seen in the Na2·E1·(ADP·Pi) simulation (Figure 6A).
Under physiological conditions, the challenge faced by the Na+/K+-pump is to effectively pick out the correct ion species from a solution much more concentrated with other types of ions. Even when the binding affinity of the ion species being transported is a few orders of magnitude higher than that of the other ion ( = 1 to 3 kcal/mol), the advantage in binding is undermined by the higher concentration of the competing ion. The pump overcomes this problem by raising the free energy barrier for the occlusion step in the presence of incorrect ions in the binding pocket, thus preventing the faulty transport (Figure 5 and Table 5). This self-correcting mechanism makes sure that the intracellular and the extracellular gates do not close (i.e., occlude), which is required for the pump cycle to go forward, unless the correct ion configuration is present at the binding pocket. Therefore, this is an inherent part of the ping-pong mechanism the pump uses to transport ions with high selectivity. Although other energetic barriers, like pathways the ions use to travel to the binding site, could contribute to the ion binding specificity, it is not likely the case here as both Na+ and K+ can access the binding pocket in the simulations of the E1 and P-E2 states without ions bound.
The concept of a self-correcting pumping mechanism sheds new light on a number of puzzling functional observations. For instance, it explains why E1 displays a strong apparent affinity for K+ in competitive binding assays (Schneeberger et al., 2001), even though the pump still binds and occludes Na+ from the cytoplasmic side. The calculations show that all the sites in the Na3·E1·(ADP·Pi) state are selective for Na+, likely because this is a partially occluded state. A wrong combination of binding site ions would not have been occluded and reached such a state. In a fully inward-open conformation, the selectivity of the sites ought to be fairly weak. The high affinity for K+ observed in E1 in these experiments is due to the backward occlusion leading to the K+-bound state E2(K2). The proposed self-correcting mechanism parallels the suggestion based on experiments that the phosphorylation and occlusion of E1 requires 3 Na+ bound, and this increases the apparent affinity for Na+ in the normal pump cycle (Schneeberger et al., 2001). The self-correcting mechanism can also account for the different effects of Na+ and K+ congeners on the release rate of the occluded 86Rb+ (a K+ congener) to the extracellular side upon the “backdoor” phosphorylation in the presence of Pi reported by Forbush (Forbush, 1987). This backward deocclusion is accelerated in the presence of Na+, and it follows a single phase, while a much slower deocclusion process with a second slow phase is observed in the presence of K+ or Rb+. According to the self-correcting mechanism, it is likely that upon the replacement of one of the two occluded 42K+ or 86Rb+ by cold K+ at one site, the pump succeeded to occlude again, resulting in the second slow phase. These explanations are sensible, in turn lending support to the proposed self-correcting mechanism. Given that the selectivity in the outward/inward facing states is after all not as strong as previously thought, such a mechanism must be in place.
This analysis suggests that an important component of the overall selectivity emerges from the increased free energy barrier associated with the occlusion process while the wrong type of ion is bound to a site. Although the pump would be kinetically more efficient by preventing the wrong ions from reaching the binding sites to begin with, the free energy calculations do not support the notion of highly selective binding sites for the open-access states. While such mechanism may seem inefficient because the pump must try to figure out that there is an issue with selectivity only after binding of several ions of the wrong type, it is important to realize that the Na+/K+-pump is not a particular fast molecular machine. The estimated turnover rate is less than 100 per second and decreases even further as the transmembrane voltage becomes more negative (Heyse et al., 1994). The implication is that the system has plenty of time, at the molecular level, to function near thermodynamic equilibrium. In other words, the pump has not been evolutionarily optimized to be a particularly fast molecular machine, but an energetically efficient one. Thus, even though the idea of a self-correcting ion selectivity mechanism seems counterintuitive and inefficient, it is consistent with the physical conditions under which the pump has to operate.
In summary, the present study highlights the tight coupling between the Na+/K+ selectivity of the binding sites, the protonation state of the coordinating residues, and the conformational state of the pump. The important functional consequence of such tight coupling is the necessity to have the correct type and number of ions in the binding pocket for the pumping cycle to proceed forward toward the occlusion step. Because the ion binding selectivity is strongly dependent upon the protein conformation, a self-correcting mechanism counteracting the effect of the ion concentration in the environment ensures an efficient function of the pump.
The crystal structures, 3WGV (Kanai et al., 2013) and 2ZXE (Shinoda et al., 2009), representing the Na+/K+-pump in its Na3·E1(ADP·Pi) and E2(K2) states were used to build the simulation systems. The structure 3WGV contains two copies of the pump assembly in the asymmetric unit. Since the structural variations between the two assemblies are minimal, only the copy including chains A, B, and G was kept. Several small molecule ligands were co-crystalized with the pump in 3WGV, including an ADP, a ion, an oligomycin A, two Mg2+, three cholesterol, four Na+, and five 1,2-diacyl-sn-glycero-3-phosphocholine molecules for each copy of the pump structure. Among these ligands, oligomycin A was not included in the simulation system. The ion was replaced by a and POPC molecules were used in place of the 1,2-diacyl-sn-glycero-3-phosphocholine. The structure of 2ZXE also contains small molecule ligands including cholesterol, K+, Mg2+, and . Similarly, a was used to replace the ion. Since the Na+/K+-pump crystal structures were obtained from different organisms, the residue numbers differ slightly. For the sake of convenience, the numbering scheme in the newly resolved 3WGV from pig kidney was adopted in this study.
After removing irregularities from the PDB files, the Na+/K+-pump subunits were capped with acetylated N-termini and amidated C-termini. The ectodomain of the β-subunit was not included in simulation systems to reduce the computational cost. The orientation was chosen to be the same as in the OPM database (Lomize et al., 2006). At this stage, the protonation states of the binding site residues were assigned with the PATCH command in CHARMM (Brooks et al., 2009). Eight different protonation states were considered in Na3·E1(ADP·Pi) state systems and both the protonated and charged D926 were included in the E2(K2) state systems. This resulted in a total of ten systems before proceeding to the next step. Table 3 shows the system names and the associated binding site protonation. We used the membrane builder module in CHARMM-GUI (Jo et al., 2007, 2008, 2009; Wu et al., 2014) to generate POPC bilayers around the pump structures. Once this was completed, each protein-membrane complex was solvated by an equal molar mixture of KCl and NaCl. The final system had a combined cation concentration of 0.3 M (i.e., [K+] = [Na+] = 0.15 M). At the end the E1 system was 84 × 110 × 158 Å3 in size and contained ~138,000 atoms, while the dimension of the E2 system was 85 × 109 × 105 Å3 with ~134,000 atoms. Each completed system was subjected to a 675-ps equilibration with reducing restraints on the heavy atoms to relax the initially uncorrelated components, followed by a 10-ns unrestrained pre-production using the simulation package NAMD2.9 (Phillips et al., 2005). After the systems were well equilibrated, they were simulated longer using the special-purpose supercomputer Anton (Shaw et al., 2009), which is designed for long time scale MD simulations.
Experimentally asparagine and glutamine are used as surrogates for protonated aspartate to study the effect of protonation. The effects of these charge-neutralizing mutants on the selectivity of the pump can be investigated computationally. The most interesting mutations in the context of this study are the single D to N mutations at the binding pocket, including residues D804, D808, and D926. The mutations were made on the crystal structures by replacing the proton on the OD2 atom in the aspartate with an –NH2 amine group. D804N and D808N mutations were also made in the outward facing P-E2 structural model (see below) to study how they affect the external K+ binding. The protonation states of the other titratable residues were determined with PROPKA. Each of these systems was equilibrated without any restraints for 40 ns. The system snapshot with the least structural deviation of the pump to the averaged structure during the simulation was used to generate the reduced system. The mutant systems are shown in Table 3.
The absolute free energy of an ion i binding to a binding site inside a protein has the following form,
The difference in the first term, , represents the nonbonded interaction component of the free energy change upon moving the ion from the bulk solution to the binding site. The subtraction in the second term, , reflects the lost of translational freedom. The translational freedom factor Ft in bulk solution can be evaluated analytically under a rigid rotor approximation (Deng et al., 2006). Its final form depends on the force constants and the equilibrium values in the distance and angle restraints applied on the ion and the surrounding protein atoms. Based on Equation (1), the selectivity of a binding site can be expressed as the binding free energy difference of two ion species. For example, the binding free energy change upon changing ion i to j at the binding site is
A negative value of indicates that ion j binds more favorably than ion i and the site is j selective. There are three terms to be evaluated in Equation (2). and are calculated using the reduce binding site model, while is computed using a water sphere with the impact from the bulk solution factored in with a boundary potential (see below).
To generate the reduced binding site, the all-atom system prepared according to the procedures above was divided into an inner region and an outer region. The inner region was defined as residues and water molecules within 15 Å to the center of mass of the bound ions. Everything within this region was treated explicitly. An extended inner region was specified by a 3-Å thick shell continuing from the inner region outwards to create a smooth spherical dielectric cavity. Water molecules in this region were removed and their impact on the binding site was included implicitly. The coordinates of atoms in these extended region and those linked to them within three atomic bonds in the inner region were held fixed during the simulations. The rest of the system was considered the outer region, in which the atoms were removed and their impact on the inner region atoms was represented by the General Solvent Boundary Potential (GSBP) in the form of a solvent-shielded static field and a solvent-induced reaction field (Im et al., 2001). The reaction field arising from changes in charge distribution in the inner region was expressed in terms of a generalized multipole expansion using 11 spherical harmonic functions. Both the solvent-shielded static field and the reaction field matrix were independent of the inner region configuration, and therefore were calculated only once before further simulations using the finite-difference Poisson−Boltzmann (PB) method with the PBEQ module (Im et al., 1998) in CHARMM. In these calculations a dielectric constant of 1 was used for the inside of the protein within the inner and outer regions, whereas the rest of the system had a dielectric constant of 80. The atomic Born radii used in the PB calculations were determined by free energy calculations in explicit solvent (Nina et al., 1997). All these reduced systems were further equilibrated for 200 ps using Langevin dynamics at 303.15 K with a friction coefficient of 5 ps−1 assigned to all non-hydrogen atoms. All bonds involving hydrogen atoms were fixed with the SHAKE algorithm (Ryckaert et al., 1977). Nonbonded interactions within 14.5 Å were accounted for explicitly, while everything else beyond this distance was treated with an extended electrostatics method (Stote et al., 1991). The simulation program CHARMM (Brooks et al., 2009) was used for the equilibration. Table 3 lists all the reduced systems.
It is known that the binding free energy calculated from FEP/MD simulations suffers from convergence issues because the residue sidechains at the binding site only sample a few of the several accessible rotameric states. This problem can be alleviated by introducing a replica-exchange scheme aiming at enhancing the sampling of sidechains (Jiang et al., 2009, 2010). This scheme allows exchange between the neighboring λ windows. Each λ window is coupled with a set of replica systems with a boosting potential of increasing strength used to accelerate the inter-conversion between sidechain rotameric states. We employed this hybrid FEP/H-REMD scheme implemented in the REPDSTR module in CHARMM (Jiang et al., 2010) to calculate the .
The boosting potential for the χ1 sidechain torsion of a binding site residue was obtained by fitting the potential of mean force as a function of the torsion χ, , with a cosine Fourier series in the form of
The angle χ1 is dihedral formed by the bonded atoms N, CA, CB, and CG. The total number of the cosine terms, N, varies from 3 to 6, depending on which one produce a better fit to the . A residue is considered in the binding site if any of its sidechain heavy atoms is within 4.5 Å to the bound Na+ or K+. Table 6 lists all the binding site residues included in the boosting potential calculations.
We performed umbrella sampling simulations to obtain the for each binding site residue. First, the entire transmembrane helix containing the residue of interest was taken from the crystal structure with its orientation kept the same as in the OPM database. This was to provide a proper secondary structure environment. Next, we replaced the sidechains of all the other residues on the helix to hydrogen atoms to remove their steric effects. An implicit planar membrane model, EEF1/IMM1 (Lazaridis, 2003; Schneeberger et al., 1999) was used in place of a membrane bilayer to solvate the helix. 72 umbrella windows along χ were set up with a window spacing of of 5°. A harmonic bias potential was applied in each window with a force constant of 100 kcal/mol/rad2. After the windows were generated, each of them was simulated for 50 ns at 303.15 K using Langevin dynamics with a 2 fs time-step. All bond lengths involving hydrogen were fixed with the SHAKE algorithm (Ryckaert et al., 1977). The cutoff distance for nonbonded interactions was set to 11 Å. The simulation program CHARMM (Brooks et al., 2009) was used to conduct the simulations. The resulting distributions of χ were unbiased using the weighted histogram analysis method (WHAM) (Kumar et al., 1992). Figure 7 shows the and the fitted curves. The fitting parameters and are given in Table 6. An appropriate boosting potential can be easily applied using the CONS DIHE command in CHARMM, with the sign in front of reversed. This will effectively cancel out the potential barrier between the rotamers of a given residue sidechain.
Before computing and , restraints were set up to restrict the translational freedom of the bound ion of interest. First, three points within the protein region were picked out. They were used along with the position of the bound ion to set up the restraints. We employed the same protocol as in the Ligand Binder module (Jo et al., 2013) in CHARMM-GUI to select the three protein anchor points, p1, p2, and p3. The relative position of the bound ions to the protein can be defined by the distance r from l1 to p1, the angle θ between l1, p1 and p2, and the torsion ψ along l1, p1, p2, and p3. The translational restraint potential took the form of
The force constant regarding distance (i.e., kr) was set to 10 kcal/mol/Å2, and the rest of the force constants were 200 kcal/mol/rad2. The equilibrium values in utrans were taken from the 200 ps equilibration of the reduced binding site.
was computed with a set of FEP/H-REMD simulations. The translational restraint, utrans, was applied to restrict the translational freedom of the ion. The Na+ ion was changed alchemically to a K+ with an alchemical coupling factor, λ. 16 λ windows were used. Each λ window included 8 replicas with the strength of boosting potential scaled from 0 to 1. Exchange attempts were made every 0.2 ps and were only allowed between neighboring replicas with different boosting potential strengths in the same λ window and between the neighboring λ windows with zero boosting potential. A total of 128 replica systems were included in the calculations. Each replica system was simulated for 2 ns. To compute at a given binding site, two sets of simulations were set up, one with Na+ at the binding site and the other with K+ for the calculation of and . The translational and orientational restraints were decoupled gradually (Beglov et al., 1994) using a coupling factor λ. The λ window and the boosting potential setup were similar to those used in the calculations. All the FEP/H-REMD calculations were performed using the REPDSTR module in CHARMM (Jiang et al., 2010). The output energies from the zero boosting potential systems were collected and processed using WHAM (Kumar et al., 1992). To compute the standard error of , we divided the trajectories into 10 blocks. The standard deviations of the averaged from all the blocks are computed and reported in Tables 4 and 5 and Figures 3–5. Cautions should be taken when comparing the calculated and experimentally measured , because the latter usually contains contributions from multiple states of the pump while the calculated is done using a single state.
The bulk system was generated by building a water sphere of 10 Å radius and centering it at the origin. The water sphere was made form pre-equilibrated water boxes with TIP3P water molecules. A Na+ ion was placed at the center of the sphere. Any water molecules overlapping with the ion were deleted. The influence of the remaining bulk was approximated by the spherical solvent boundary potential (Beglov et al., 1994). The system was equilibrated for 200 ps at 303.15 K. Other simulation options were kept the same as described in the reduced binding site system. During the equilibration the position of the ion was restrained using a weak harmonic bias potential with a force constant of 0.5 kcal/mol/Å2. Once equilibrated, the Na+ was gradually changed to a K+ using the PERT module in CHARMM with 11 λ windows. The simulation time for each λ window was 1 ns. The resulting is 18.34 kcal/mol after unbiasing the energy outputs with WHAM. The calculated for all the binding sites in the systems summarized in Table 3 are shown in Table 4.
To further confirm the protonation state of D804, we evaluated its pKa shift with additional simulations in explicit solvent using the following formula:
and are the free energy change of aspartate deprotonation in the protein environment and in bulk water, respectively. The reduced system of E1_S1 (Table 3) was used to compute the deprotonation free energy of D804 at the ion binding site (). To calculate the , an aspartic acid residue with an acetylated N-terminus and an amidated C-terminus was put into a pre-equilibrated water sphere of 15 Å radius. As in the calculations of and , the impact of the bulk solution beyond the current system was incorporated by the spherical solvent boundary potential (Beglov et al., 1994). Alchemical FEP calculations were carried out in both systems. The λ windows were evenly spaced to gradually deprotonate the aspartate. The numbers of windows used in the and calculations were 24 and 10 respectively. Each window was equilibrated for 200 ps and contained 5 ns of sampling. The calculated and are −44.8 kcal/mol and −51.9 kcal/mol, respectively. Using Equation (5), the ΔpKa is 5.1, meaning the pKa of D804 is right shifted by 5.1 pK unit when it is in the Na+/K+ pump ion binding site. The final calculated pKa of D804 is 4.1 + 5.1 = 9.2 with 4.1 being the pKa value of an aspartate in bulk water (Berg et al., 2002), reinforcing the conclusion that D804 is protonated in the E1 state of the pump.
Three conformational transition pathways were generated based on the shared homology between the Na+/K+-pump and the Ca2+ SERCA pump. One of such pathways connects states E2(K2) and P-E2·K2. The second connects states E2(K2) and K2·E1, and the third describes the transition between states Na3·E1-P and P-E2·Na3. First, A Cα-atom-only transition pathway for each of these was generated using the ANMPathway online server (http://anmpathway.lcrc.anl.gov/anmpathway.cgi) (Das et al., 2014) based on the SERCA pump crystal structures including 3B9B (Olesen et al., 2007), 1WPG (Toyoshima et al., 2004), and 1VFP (Toyoshima et al., 2004). Since both the Na+/K+-pump and the Ca2+ SERCA pump are P-type ATPases, it is reasonable to assume that they have similar pumping mechanisms and they share the same set of states along the pump cycle. Although crystal structures of the Na+/K+-pump are scarce, multiple high-resolution structures of the SERCA pump in different states are available (Karlsen et al., 2016). Among them, 3B9B represents the outward facing P-E2 state. All the other SERCA pump structures were aligned based on their transmembrane region Cα positions to the Na+/K+-pump and the Cα RMSD were computed to find those that resemble the structural states captured by the Na+/K+-pump crystal structures 3WGV (Kanai et al., 2013) and 2ZXE (Shinoda et al., 2009). The two SERCA pump structures representing states Na3·E1ATP (3WGV) and E2(K2) (2ZXE) are 1VFP and 1WPG, respectively. Each generated coarse-grained pathway is made of a sequence of structural snapshots containing Cα atoms only. These snapshots are called images and are distributed at equal RMSD intervals. The images from the coarse-grained pathways were then used in the all-atom targeted molecular dynamics (TMD) simulations to guide each transition. The starting system in the first and second pathway TMD simulations was taken from the well-equilibrated MD simulation system E2_S1 starting from the crystal structure 2ZXE. The third pathway was realized using the MD simulation system E1_S1 built from the crystal structure 3WGV, with the catalytic D369 phosphorylated. The protocol of the TMD simulations followed those published before (Castillo et al., 2015)
Binding site ion selectivity was evaluated for nine systems representing different stages along the pump cycle. Among these systems are well-defined states: P-E2·K2, E2(K2), K2·E1, Na3·E1(ADP·Pi), and P-E2·Na3. The structures used for them are ① the generated P-E2·K2 model, ③ the equilibrated crystal structure 2ZXE, the equilibrated crystal structure 3WGV with ⑦ bound Na+ and with ⑤ bound Na+ replaced by K+ at sites I and II, and ⑨ the generated P-E2·Na3 model. Several intermediate states were also included. Intermediate state ② is between states ① P-E2·K2 and ③ E2(K2), taken from the TMD simulations of the first path after image 35. State ④ was taken at image 33 along the transition from state ③ E2(K2) to ⑤ K2·E1. Along the same vein, state ⑧ represented the intermediate at image 29 between states ⑦ Na3·E1(ADP.Pi) and ⑨ P-E2·Na3. The intermediates were taken near the midpoints of the transitions. An additional state ⑥ K2·E1* with altered binding site residue protonation was also included to reflect the protonation state change upon the E2 to E1 transition. All the states included in the selectivity calculations and their relative placement along the pump cycle can be seen in Figure 5.
To evaluate the selectivity at a given binding site in a system, three simulations were performed with slightly varied Lennard-Jones (LJ) parameters for the ion of interest. In the first simulation, the parameters of the ion were unaltered. A linear interpolation of the LJ and the NBFIX parameters between Na+ and K+ was made. In the second simulation, the LJ and NBFIX parameters of the ion were replaced by those of a Na/K hybrid at the middle point of the interpolation. In the last simulation, the ion’s parameters were completely changed to represent the other ion type, either K+ or Na+, depending on the starting ion type. 10 ns trajectories were generated for each simulation using the simulation package OpenMM6.2 (Eastman et al., 2013). Energies were calculated from each trajectory using two parameter sets. The difference was fed into the WHAM equation (Kumar et al., 1992) and the free energy changes upon mutating the ion to the intermediate as well as upon mutating the intermediate to the other ion type were solved self-consistently. The sum of the two free energies gives the . The same strategy was used to compute the in water ( kcal/mol). The binding free energy difference is given by . The selectivity of binding site with different occupancies was also included in the calculations. The are shown in Table 5.
The initial relaxation and the restraint-free equilibration of the membrane pump systems were performed using the NAMD2.9 (Phillips et al., 2005) simulation package with the input scripts from the CHARMM-GUI membrane builder module. The NAMD simulation temperature was set to 303.15 K using Langevin Dynamics with a damping coefficient of 10 ps−1 during the relaxation and 1 ps−1 during the restraint-free equilibration. The van der Waals interactions were smoothly switched off at 10–12 Å by a force switching function (Steinbach et al., 1994) and the electrostatic interactions were calculated using the particle-mesh Ewald method with a mesh size of ~1 Å.
After a short equilibration with NAMD, each all-atom system listed in Table 3 were subjected to a hundred ns scale long production without any restraints using the special-purpose supercomputer Anton (Shaw et al., 2009). The volume of the periodic cell was kept constant and the temperature was set to 303.15K using the Nosé-Hoover thermostats (Martyna et al., 1992). The lengths of all bonds involving hydrogen atoms were constrained using M-SHAKE (Kräutler et al., 2001). The cutoff of the van der Waals and short-range electrostatic interactions was set to an optimal value suggested by the Anton guesser script, guess_chem. Long-range electrostatic interactions were evaluated using the k-space Gaussian split Ewald method (Shan et al., 2005) with a 64 × 64 × 64 mesh. The integration time step was 2 fs. The r-RESPA integration method (Tuckerman et al., 1992) was employed and long-range electrostatics were evaluated every 6 fs.
Simulations used to compute ion selectivity along the pump cycle were conducted using the simulation package OpenMM6.2 (Eastman et al., 2013). The constant pressure and temperature (NPT) ensemble was used for these simulations. The pressure was maintained at 1 atm with a Monte Carlo barostat, which attempts to adjust the system volume every 0.2 ps. The Langevin dynamics algorithm with a 1.0 ps−1 friction coefficient was used to hold the simulation temperature, which is at 303.15 K. The lengths of all bonds involving hydrogen were fixed and the integration time step was set to 2 fs. A force switching function was applied from 10 to 12 Å to gradually turn off the van der Waals interactions and the particle-mesh Ewald method with an error tolerance of 0.0005 was used to evaluate the electrostatic interactions.
The same force field parameters were used in all other simulations in this study. The PARAM27 all-atom force field of CHARMM (MacKerell et al., 1998) with a modified version of dihedral cross-term correction (Mackerell et al., 2004) was used for the protein and the C36 lipid force field (Klauda et al., 2010) was used for POPC. Water molecules were modeled with the TIP3P potential (Jorgensen et al., 1983).
Biochemical aspects of active transportAnnual Review of Biochemistry 36:727–756.https://doi.org/10.1146/annurev.bi.36.070167.003455
Evolution of substrate specificities in the P-type ATPase superfamilyJournal of Molecular Evolution 46:84–101.https://doi.org/10.1007/PL00006286
Finite representation of an infinite bulk system: Solvent boundary potential for computer simulationsThe Journal of Chemical Physics 100:9050–9063.https://doi.org/10.1063/1.466711
Biochemistry (5th ed)W.H. Freeman.
Mechanism of potassium ion uptake by the Na(+)/K(+)-ATPaseNature Communications 6:7622.https://doi.org/10.1038/ncomms8622
The Pymol Molecular Graphics SystemSchrödinger, Palo Alto, CA.
Calculation of standard binding free energies: Aromatic molecules in the T4 Lysozyme L99A mutantJournal of Chemical Theory and Computation 2:1255–1273.https://doi.org/10.1021/ct060037v
OpenMM 4: A reusable, extensible, hardware independent library for high performance molecular simulationJournal of Chemical Theory and Computation 9:461–469.https://doi.org/10.1021/ct300857j
The rapid-onset dystonia parkinsonism mutation D923N of the Na+, K+-ATPase alpha3 isoform disrupts Na+ interaction at the third Na+ siteJournal of Biological Chemistry 285:26245–26254.https://doi.org/10.1074/jbc.M110.123976
Rapid release of 42K or 86Rb from two distinct transport sites on the Na,K-pump in the presence of Pi or vanadateThe Journal of Biological Chemistry 262:11116–11127.
Isolation, crystallizationand crystal structure determination of bovine kidney na(+),k(+)-atpaseActa Crystallographica Section F Structural Biology Communications 72:282–287.https://doi.org/10.1107/S2053230X1600279X
Partial reactions of the Na,K-ATPase: determination of rate constantsThe Journal of General Physiology 104:197–240.https://doi.org/10.1085/jgp.104.2.197
Continuum solvation model: Computation of electrostatic forces from numerical solutions to the Poisson-Boltzmann equationComputer Physics Communications 111:59–75.https://doi.org/10.1016/S0010-4655(98)00016-2
Generalized solvent boundary potential for computer simulationsThe Journal of Chemical Physics 114:2924–2937.https://doi.org/10.1063/1.1336570
Modulation of the Na,K-pump function by beta subunit isoformsThe Journal of General Physiology 103:605–623.https://doi.org/10.1085/jgp.103.4.605
Computation of absolute hydration and binding free energy with free energy perturbation distributed replica-exchange molecular dynamics (FEP/REMD)Journal of Chemical Theory and Computation 5:2583–2588.https://doi.org/10.1021/ct900223z
Free energy perturbation hamiltonian replica-exchange molecular dynamics (fep/h-remd) for absolute ligand binding free energy calculationsJournal of Chemical Theory and Computation 6:2559–2565.https://doi.org/10.1021/ct1001768
CHARMM-GUI Ligand Binder for absolute binding free energy calculations and its applicationJournal of Chemical Information and Modeling 53:267–277.https://doi.org/10.1021/ci300505n
Software news and updates - charnim-gui: A web-based grraphical user interface for charmmJournal of Computational Chemistry 29:1859–1865.https://doi.org/10.1002/jcc.20945
Structure–function relationships of Na+, K+, ATP, or Mg2+ binding and energy transduction in Na,K-ATPaseBiochimica Et Biophysica Acta (BBA) - Bioenergetics 1505:57–74.https://doi.org/10.1016/S0005-2728(00)00277-2
Comparison of simple potential functions for simulating liquid waterThe Journal of Chemical Physics 79:926–935.https://doi.org/10.1063/1.445869
How to compare, analyze, and morph between crystal structures of different conformations: The P-Type ATPase exampleMethods in Molecular Biology 1377:523–539.https://doi.org/10.1007/978-1-4939-3179-8_43
Update of the CHARMM all-atom additive force field for lipids: validation on six lipid typesThe Journal of Physical Chemistry B 114:7830–7843.https://doi.org/10.1021/jp101759q
THE weighted histogram analysis method for free-energy calculations on biomolecules. I. The methodJournal of Computational Chemistry 13:1011–1021.https://doi.org/10.1002/jcc.540130812
Asp804 and Asp808 in the transmembrane domain of the Na,K-ATPase alpha subunit are cation coordinating residuesJournal of Biological Chemistry 271:29682–29687.https://doi.org/10.1074/jbc.271.47.29682
Glutamic acid 327 in the sheep alpha 1 isoform of Na+,K(+)-ATPase stabilizes a K(+)-induced conformational changeThe Journal of Biological Chemistry 270:2993–3000.https://doi.org/10.1074/jbc.270.7.2993
Effective energy function for proteins in solutionProteins: Structure, Function, and Genetics 35:133–152.https://doi.org/10.1002/(SICI)1097-0134(19990501)35:2<133::AID-PROT1>3.0.CO;2-N
Effective energy function for proteins in lipid membranesProteins: Structure, Function, and Genetics 52:176–192.https://doi.org/10.1002/prot.10410
OPM: orientations of proteins in membranes databaseBioinformatics 22:623–625.https://doi.org/10.1093/bioinformatics/btk023
All-atom empirical potential for molecular modeling and dynamics studies of proteinsThe Journal of Physical Chemistry. B 102:3586–3616.https://doi.org/10.1021/jp973084f
K+ congeners that do not compromise Na+ activation of the Na+,K+-ATPase: hydration of the ion binding cavity likely controls ion selectivityThe Journal of Biological Chemistry 290:3720–3731.https://doi.org/10.1074/jbc.M114.577486
Nosé–Hoover chains: The canonical ensemble via continuous dynamicsThe Journal of Chemical Physics 97:2635–2643.https://doi.org/10.1063/1.463940
Energy costs of ion pumping by animal tissuesThe Journal of Nutrition 115:1374–1382.
Sodium and proton effects on inward proton transport through Na/K pumpsBiophysical Journal 106:2555–2565.https://doi.org/10.1016/j.bpj.2014.04.053
Atomic radii for continuum electrostatics calculations based on molecular dynamics free energy simulationsThe Journal of Physical Chemistry B 101:5239–5248.https://doi.org/10.1021/jp970736r
PROPKA3: Consistent treatment of internal and surface residues in empirical pKa predictionsJournal of Chemical Theory and Computation 7:525–537.https://doi.org/10.1021/ct100578z
Flexibility of an active center in sodium-plus-potassium adenosine triphosphataseThe Journal of General Physiology 54:306–326.https://doi.org/10.1085/jgp.54.1.306
Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanesJournal of Computational Physics 23:327–341.https://doi.org/10.1016/0021-9991(77)90098-5
Ion selectivity of the cytoplasmic binding sites of the Na,K-ATPase: II. Competition of various cationsJournal of Membrane Biology 179:263–273.https://doi.org/10.1007/s002320010051
Gaussian split Ewald: A fast Ewald mesh method for molecular simulationThe Journal of Chemical Physics 122:054101.https://doi.org/10.1063/1.1839571
Millisecond-scale molecular dynamics simulations on AntonProceedings of the Conference on High Performance Computing Networking.https://doi.org/10.1145/1654059.1654126
Effects of ATP and protons on the Na : K selectivity of the (Na+ + K+)-ATPase studied by ligand effects on intrinsic and extrinsic fluorescenceBiochimica Et Biophysica Acta (BBA) - Biomembranes 601:386–402.https://doi.org/10.1016/0005-2736(80)90543-X
New spherical-cutoff methods for long-range forces in macromolecular simulationJournal of Computational Chemistry 15:667–683.https://doi.org/10.1002/jcc.540150702
On the treatment of electrostatic interactions in biomolecular simulationJournal De Chimie Physique Et De Physico-Chimie Biologique 88:2419–2433.https://doi.org/10.1063/1.41345
Effect of extracellular pH on presteady-state and steady-state current mediated by the Na+/K+ pumpJournal of Membrane Biology 198:65–76.https://doi.org/10.1007/s00232-004-0660-4
CHARMM-GUI Membrane Builder toward realistic biological membrane simulationsJournal of Computational Chemistry 35:1997–2004.https://doi.org/10.1002/jcc.23702
Protonation of key acidic residues is critical for the K⁺-selectivity of the Na/K pumpNature Structural & Molecular Biology 18:1159–1163.https://doi.org/10.1038/nsmb.2113
Nir Ben-TalReviewing Editor; Tel Aviv University, Israel
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Binding site protonation and self-correcting occlusion control the Na+/K+-pump selectivity" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Erik Lindahl (Reviewer #3) and Nir Ben-Tal (Reviewer #1), who is a member of our Board of Reviewing Editors, and the evaluation has been overseen by Richard Aldrich as the Senior Editor.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
The study examines the physicochemical determinants of binding site selectivity in the Na+/K+-ATPase using molecular dynamics simulations. Various conformation of the ATPase, corresponding to different states along the pump cycle are examined, including the crystal structures of its two main conformations, as well as several models of intermediate states based on the Ca2+ SERCA pump, and models derived from interpolating between well-defined states. Various protonation states of titratable residues in the pump's binding sites are energetically examined using molecular dynamic simulations as different resolutions. The results suggest that the protonation states of three aspartate residues in the binding sites (D804, D808 and D926) alter upon switching between conformations. The results indicate that the different selectivity of the binding sites emerges from the change in the protonation states. Using reduced structural models of the binding sites, this hypothesis is supported by calculating the difference in free energy between the pump bound to Na+ vs. K+ ions along the transport cycle. The authors conclude that initially the binding sites are only weakly selective for either Na+ or K+, and selectivity emerges as collective effect due to the (indirect) interaction of new ions with the ones that are already bound (a knock-on effect); the pump shifts towards K+ selectivity in the outward open state, and towards Na+ selectivity in the inward open state. Binding of the wrong type of ion increases the free energy barrier towards occlusion, preventing the transport cycle to proceed, thereby preventing the pump from operating in the opposite (wrong) direction.
The manuscript addresses an interesting and important question. It proposes a simple mechanism that dictates the shift in the selective properties of the Na+/K+-pump binding sites. However, based on the current draft of the manuscript it is difficult to decide whether the work is novel and solid enough to justify publication in eLife. The summary reflects a long discussion among the reviewers in an effort to resolve this. The author should submit a revised draft only if they are certain that they can address all the issues raised.
1) It is difficult to figure out which structures the various steps in Figure 5 are based on. Perhaps:
State 1 and 9 correspond to 3B9B (open SERCA);
State 2 calculated transition to;
State 3 corresponds to 1WPG (occluded SERCA) – but why not 2ZXE (occluded NKA)?
State 4 calculated transition to;
State 5, unclear: called K2E1. Maybe 1VFP with different protonation pattern and K instead of Na?
State 6 corresponds to 1VFP (Ca occluded SERCA) – but with K rather than Na?
State 7 corresponds to 3WGV (occluded NKA);
State 8 calculated transition to 9;
The authors should make perfectly clear which structures are used for which states.
2) The pump cycle should be based on all the relevant structures. Structures have been published of SERCA in E2 (3W5C) and the open ready-to-bind Ca (4H1W) in 2013. These important steps are missing here, or guesses are attempted based on the Ca-occluded 1VFP structure, but that did not predict the structure for H and Ca exchange. The structures are not mentioned. (A nice and recent overview of P-type structures can be found here http://www.ncbi.nlm.nih.gov/pubmed/26695058).
Also, it would be preferable to use 4HYT (ouabain-bound NKA, mentioned as recent but it is from 2013) rather than 3B9B. At least the authors could show an overlay of the ion binding sites from 4HYT and their model.
3) Standard errors estimates should be added to the energy values, and conclusions that have been based on energy values with too large 'error bars' should be reconsidered.
4) If the calculations correctly describe the binding free energy of the pump cycle, they would also need to explain the ion release on the opposite side. Can they show this?
5) The simulations and free energy calculations assume a starting point in which all the pump's binding sites are occupied by a certain set of ions. The authors explicitly claim that selectivity emerges from the increased free energy barrier associate with binding of the wrong ion to a site. However, from a biological point of view to argue that selectivity is governed solely by such mechanism is puzzling, for it is inefficient. The mechanism proposed here beautifully explains the different selective properties of the pump's binding sites in its two conformations, and it is possibly a "last resort mechanism", designed, as the authors claim, to prevent the transport of the wrong kind of ion if it reaches the binding site. However one would expect that for the most part, under physiological conditions, the pump will prevent the wrong ions from reaching the binding site to begin with. The anticipation is that a fast machine, such as this pump, should follow a particularly efficient mechanism, rather than the one proposed here, which figures out that there is issue with selectivity only after binding of several ions of the wrong type. The manuscript does not account for the ions transition into the binding sites and other possible energetic barriers they might encounter on the way. This issue should, at the very least, be discussed.
6) In view of item 5, the authors should revise the title to reflect the fact that selectivity is presumably determined primarily by other mechanisms, and the suggested mechanism here is only a 'safe guard'.https://doi.org/10.7554/eLife.16616.016
- Huan Rui
- Benoît Roux
- Pablo Artigas
- Benoît Roux
- Pablo Artigas
- Benoît Roux
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
This work was supported by grant U54-GM087519 (to BR) and R15(NS081570-01A to PA) from NIH, and by grant MCB-1515434 (to PA and BR) from NSF. Anton computer time was provided by the National Center for Multiscale Modeling of Biological Systems (MMBioS) through Grant P41GM103712-S1 from the National Institutes of Health and the Pittsburgh Supercomputing Center (PSC). The Anton machine at PSC was generously made available by D.E Shaw Research.
- Nir Ben-Tal, Reviewing Editor, Tel Aviv University, Israel
© 2016, Rui et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.