Deprotonation by Dehydration: The Origin of Ammonium Sensing in the AmtB Channel

The AmtB channel passively allows the transport of NH4 + across the membranes of bacteria via a “gas” NH3 intermediate and is related by homology (sequentially, structurally, and functionally) to many forms of Rh protein (both erythroid and nonerythroid) found in animals and humans. New structural information on this channel has inspired computational studies aimed at clarifying various aspects of NH4 + recruitment and binding in the periplasm, as well as its deprotonation. However, precise mechanisms for these events are still unknown, and, so far, explanations for subsequent NH3 translocation and reprotonation at the cytoplasmic end of the channel have not been rigorously addressed. We employ molecular dynamics simulations and free energy methods on a full AmtB trimer system in membrane and bathed in electrolyte. Combining the potential of mean force for NH4 +/NH3 translocation with data from thermodynamic integration calculations allows us to find the apparent pKa of NH4 + as a function of the transport axis. Our calculations reveal the specific sites at which its deprotonation (at the periplasmic end) and reprotonation (at the cytoplasmic end) occurs. Contrary to most hypotheses, which ascribe a proton-accepting role to various periplasmic or luminal residues of the channel, our results suggest that the most plausible proton donor/acceptor at either of these sites is water. Free-energetic analysis not only verifies crystallographically determined binding sites for NH4 + and NH3 along the transport axis, but also reveals a previously undetermined binding site for NH4 + at the cytoplasmic end of the channel. Analysis of dynamics and the free energies of all possible loading states for NH3 inside the channel also reveal that hydrophobic pressure and the free-energetic profile provided by the pore lumen drives this species toward the cytoplasm for protonation just before reaching the newly discovered site.


Introduction
The transport of (NH 4 þ ) ammonium and/or (NH 3 ) ammonia (we will refer to both of these species together as Am) across biological membranes is a homeostatic necessity in both prokaryotes and eukaryotes [1]. In the case of many different plants, bacteria, and fungi, Am serves as a readily available nitrogen source for biosynthetic purposes. On the other hand, at high concentrations, it becomes cytotoxic, especially in animal cells. The family of Am transport proteins-ammonium transporters (Amt) in plants and bacteria, methylamine permeases (MEP) in yeast, and rhesus (Rh) proteins in animals-serves to facilitate the permeation of Am across the membrane. Plant [2][3][4][5] and yeast [6,7] Amt/ MEPs as well as many bacterial [8][9][10] Amts take in Am in a membrane electrochemical potential-dependent manner in order to utilize it. In humans, the related Rh proteins are split into two groups: erythroid (RhAG, RhD, and RhCE)-expressed on the erythrocyte surface [11,12] where they perform immunogenic and structural roles, and nonerythroid (RhCG, RhBG, and RhGK)-expressed in the kidneys, liver, and testes where they aid in disposal of ammonium and regulation of pH [13,14].
Many years of study have shown that while members of the Am transporter family share homologous sequences and structures, it does not necessarily follow that they conduct Am using the same mechanism [15]. Whether particular members of the family transport Am in its ionic (NH 4 þ ) or ''gas'' (NH 3 ) form remains a subject of debate [15]. It is also suggested that some members of the family transport H þ and NH 3 in a coupled fashion [4,10,16]. Recently, X-ray diffraction studies [17,18] have revealed the atomic structure of the bacterial Am transporter, AmtB, from Escherichia coli, providing a quantum leap in our understanding of the permeation of Am in the form of NH 3 . Even though its tomato plant homolog, LeAMT1;2 does not appear to share the NH 3 transport mechanism, the fact that AmtB's human (kidney) homolog does transport NH 3 rather than NH 4 þ means that these new crystallographic structures will provide an excellent starting point for understanding the structural basis of Am transport via its ''gas'' intermediate [15] as it pertains to human physiology. In fact, the newly available structural information on AmtB has already spurred comparative structural and functional modeling of human Rh proteins [19].
In the membrane, AmtB exists as a homo-trimer, with each monomer consisting of a right-handed eleven-helix bundle [17,18]. The center of each monomer forms a narrow hydrophobic pore connecting cytoplasmic and periplasmic depressions (vestibules). The diffraction studies revealed a NH 4 þ binding site (named ''Am1'' [17]) on the periplasmic depression where the cation donates a hydrogen bond to the hydroxyl oxygen of a Ser residue (S219) and is putatively stabilized by cation-p interactions with aromatic residues (W148 and F107). Though the narrow hydrophobic pore (lumen) connecting the inner and outer vestibules of the Xray structure was seen to accommodate three poorly ordered NH 3 binding sites (named Am2, Am3, and Am4, respectively [17]), the opening from the periplasmic vestibule to the pore was apparently blocked by two stacked Phe sidechains (F107 and F215), rendering the means of NH 4 þ deprotonation and translocation unclear. The fact that the pore was lined by two His residues (H168 and H318) hydrogen bonded to one another led to the suggestion that NH 4 þ can penetrate deep into the pore up to the middle site, Am3, before undergoing His-mediated deprotonation to form NH 3 [17]. However, no site or putative mechanism for reprotonation of NH 3 to form a cation was found. Nor was any binding site for NH 4 þ found in the cytoplasmic vestibule. It was also found very unlikely that water occupies the hydrophobic pore, further supporting the notion that AmtB transports NH 3 instead of NH 4 þ . Since AmtB passively transports Am along an electrochemical gradient and was not seen to undergo conformational changes upon occupying different permeant species, it was deemed safe [17] to consider the family of Amts as ''channels.' ' The new structural knowledge of AmtB has facilitated computational studies [18,[20][21][22][23] aimed at understanding aspects of Am transport through the channel. Of particular interest in these studies is an Asp residue (D160) on the periplasmic side of the channel, which is highly conserved across Amt/MEP/Rh proteins. Mutational studies have shown this residue to be functionally crucial for Amts as well as MEP and Rh proteins [24,25]. Although it was proposed that D160 likely participates in the periplasmic binding site for NH 4 þ [24], X-ray structures revealed that the carboxylate side chain does not participate in the Am1 binding site, is sequestered from water, and interacts closely with amide N-H groups of an outer loop (connecting helices M4 and M5), hence playing an indispensable structural role in luring NH 4 þ from the periplasm by forcing it to interact with outwardly directed carbonyl oxygens [17,18]. A computational free energy perturbation study by Luzhkov et al. [22] demonstrated that the charged form of D160 is stabilized by 0.3-5.1 pK a units in the absence of a cation at site Am1. In the presence of a cation at Am1, the charged form of D160 was seen to be even more stable (by 9.2 pK a units). Performing the alchemical mutation, D160N, in agreement with experimental mutational studies [24,25] gave rise to a loss of binding for NH 4 þ at site Am1, suggesting that the presence of D160 stabilizes cation binding via throughspace electrostatic interactions. However, the simulations of Luzhkov did not evaluate the structural role of D160 in facilitating cation sequestration using backbone carbonyls in the M4-M5 loop as suggested by structural studies [17,18]. Such an evaluation would, necessarily, involve long simulations of the mutant channel. Thus, it is probably most accurate to say that D160 plays both an important role in stabilizing cation binding at site Am1 as well as stabilizing the M4-M5 loop for recruiting NH 4 þ from the bulk, and that its structural and cation-stabilizing roles likely cannot be considered separately. Computational studies performed by Lin et al. and Nygaard et al. led to the suggestion of an entirely different role for D160 as a proton acceptor from NH 4 þ , with either water [20] or the carbonyl oxygen of A162 [23] as an intermediary, based on qualitative inspection of molecular dynamics (MD) trajectories. The simulation studies of Luzhkov

Results/Discussion
To answer these questions, we performed MD simulations on a full AmtB trimer embedded in a palmitoyloleoylphosphatidylcholine (POPC) bilayer ( Figure 1)

Author Summary
Selective flow of ammonium manifests itself in a unique way in the case of the ammonium channel, AmtB, allowing it to interact closely with cytoplasmic signal transduction proteins in order to ''sense'' the presence of extracellular ammonium. Although it is well known that AmtB transports ammonia (NH 3 ) rather than ammonium ion (NH 4 þ ), it is unclear from the channel's atomic structure exactly where and how, along its pathway toward the cytoplasm, NH 4 þ becomes deprotonated to form NH 3 , and reprotonated on the cytoplasmic end of the channel to form NH 4 þ to enter the cell. We use a combination of molecular dynamics simulation techniques to glean the thermodynamics associated with these key events in ammonium translocation. Our findings provide a novel perspective on how this family of channels indirectly controls ammonium protonation-by directly controlling its hydration. Such a perspective should lend new insight to interpretations of experimental data, and could possibly lead to new strategies in an envisioned future for the design of nanopores that can control the protonated state of permeant species. and Methods). Although such a large biomembrane/electrolyte system was expected to require large equilibration times [26][27][28][29] due to slow ion-membrane adsorption, it allowed for a more large-scale view of NH 4 þ recruitment at the periplasmic and cytoplasmic vestibules. In addition, this system allowed us to determine ionic densities as a function of the transport (z-) axis from which we may derive the free energetics describing Am transport all the way from the bulk phase into the channel interior.

Vestibular NH 4 þ Recruitment and Binding
Analysis of our simulated system showed that the large dipole moment (;2,000 Debye-see Figure S1A) of the AmtB trimer has an interesting impact on the local ionic concentration at the membrane/electrolyte interface (Figures 2, S1, and S2), which comprises a dielectric response to the protein's electric field. While NH 4 þ adsorption at the membrane surface was not affected significantly, Cl À adsorption was. Overall, the local concentration of cation near the periplasmic vestibule (just beneath the lipid headgroup phosphate plane) was observed to be ;16 times greater than that near the cytoplasmic vestibule. On the other hand, Cl À was very unlikely to approach either vestibule, and its equilibrium concentration was ;3 times smaller on the periplasmic side than on the cytoplasmic side ( Figure S2). Hence, without any external electrochemical gradient, we see that the channel is built for recruiting cations from the periplasm at neutral pH on a mesoscopic scale. In addition we see that the channel inhibits anionic adsorption to the periplasmic membrane surfaceanother mechanism of enhancing cation recruitment.
To glean information about cation interaction with the channel interior, as well as specific binding interactions, we combined data from local z-dependent ion density profiles with umbrella sampling data (Materials and Methods) to calculate the potential of mean force (PMF) governing the  . Partial Densities across the Simulation Cell Shows the system structure as a function of the transport axis, z (z ¼ 0 coincides with the center of mass of the trimer backbone). (A) Major system components (lipid, water, and protein). The partial density of the phosphate group is shown as a reference for roughly identifying the boundaries of the membrane. A black, vertical dotted line is drawn at the position of the peak in each leaflet's phosphate density. (B) NH 4 þ and Cl À densities. A black vertical dotted line is drawn at leaflet phosphate positions. The average concentration of NH 4 þ and Cl À in the bulk is ''measured'' to be 158 6 5 mM, based on this plot. Cl À binding is enhanced on the cytoplasmic (z , 0) surface and inhibited on the periplasmic (z . 0) surface due to the ''external'' field created by the channel (dipole moment symbol) as seen from the Cl À density maxima at ;62.7 nm. NH 4 þ binding (see density maxima at ;61.7 nm) is not as seriously affected, although its local concentration at ;62.7 nm (around the Cl À binding positions) is enhanced near the periplasmic (compared with the cytoplasmic) membrane surface. (C) NH 4 þ binding to the membrane surface involves penetration, dehydration, and hydrogen bonding with carbonyl and ester oxygens near the glycerol backbone of lipid molecules as seen in previous work [29]. doi:10.1371/journal.pcbi.0030022.g002 transport of a single NH 4 þ molecule along the channel (z-) axis. The PMF profile, shown in Figures 3 and S3, details the free energy required to bring a single NH 4 þ molecule to any point within the channel from the bulk electrolyte on either side of the membrane. It clearly shows two free energy minima (binding sites)-one on either side of the membrane. While approaching either site from the bulk, the PMF indicates that NH 4 þ overcomes a very small barrier (;2 RT on either side) at z ' 6 2nm (z ¼ 0 nm is taken to be the center of mass of the trimer backbone). The site on the periplasmic side (at z ¼ 1.07 nm) of the membrane (labeled Am1 in Figure 3) is at precisely the location of site Am1 in the X-ray structures [17,18], although there are slight differences in the specific interactions. At this site, NH 4 þ forms hydrogen bonds with the backbone carbonyl oxygen of A162 and the sidechain hydroxyl oxygen of S219. The ''other half'' of the cation donates hydrogen bonds to a shell of water molecules, which is stabilized, in turn, by hydrogen bonds to the polar sidechain of Q104. Just below the cation, F107 provides a ''floor'' for the Am1 binding site. We also note a feature in the free-energy profile (circled in blue at z ¼ 1.342 nm) slightly above what is normally referred to as site Am1. This feature shows a slightly lower free energy than Am1, and corresponds to interaction with backbone carbonyl oxygens of D160, F161, and A162, which usher NH 4 þ into Am1. Figure 3 shows that the carboxylate of D160 accepts hydrogen bonds from the hydroxyl and amide nitrogen of T165, and from the backbone amide nitrogen from G164 and G163. These persistent interactions support D1609s structural role in orienting the carbonyl groups of the M4-M5 loop to recruit NH 4 þ from the periplasm. This is not to say that D160 does not stabilize the cation at the lower site, Am1, via through-space electrostatic interactions as suggested by Luzhkov et al. [22], but that its loop-stabilizing role (for recruitment of cations) is likely not separable from its electrostatic role. Nygaard et al. noticed in their simulations that NH 4 þ can occupy different substates (so-called Am1a or Am1b) within the Am1 binding site [23]. While, indeed, these substates may exist and contribute to the free-energy profile in Figure 3, the fact that the Am1 free-energy well is rather smooth implies that the free-energy barriers and differences between any set of substates would not be observable in a macroscopic (experimental) binding assay. In fact, the barrier (and free-energy difference) between Am1 and the feature circled in blue in Figure 3 is less than 1 RT, which implies that both features should be considered together in any discussion of ion recruitment and binding. Thus, it would make sense, from a free-energetic perspective, to refer to both of these ''substates,'' together, as the periplasmic binding site. The PMF in Figure 3 also reveals a cation binding site on the cytoplasmic side of the membrane (labeled Am5 at z ¼ À1.37 nm in Figure 3) whose existence was previously unknown according to the latest diffraction and computational studies. A NH 4 þ molecule at this new site, Am5, shows an equivalent level of stability to one bound at Am1, is mostly hydrated (usually by three to four water molecules), and donates strong hydrogen bonds to the carboxyl oxygen of D313 and the hydroxyl oxygen of S263. In units of RT ¼ 2.48 kJ/mol at 298 K. The profile was shifted such that it takes on a value of zero in the bulk. The origin is the same as in Figure 2, with periplasm (z . 0) and cytoplasm (z , 0) separated by the membrane (membrane phosphate positions are marked with vertical dotted lines). Am (de)protonation regions are marked with green vertical bars, whose widths give the uncertainty in the regions. NH 3 binding sites are marked with red vertical bars. NH 4 þ binding sites are marked Am1 (z ¼ 1.07 nm) and Am5 (z ¼À1.37 nm), along with an additional binding feature, circled in blue (z ¼ 1.32 nm). The bound state of NH 4 þ at each of the three sites is shown in molecular detail (labeled panels). All specific interactions within each site are described in the text. doi:10.1371/journal.pcbi.0030022.g003 A global view of the PMF for NH 4 þ translocation across the hydrophobic pore between the inner and outer vestibules reveals that NH 4 þ permeation would be extremely rare ( Figure S3). The free energy barrier for NH 4 þ entering the lumen from either the periplasm or cytoplasm is prohibitive (.50 RT to reach the NH 3 sites Am2, Am3, or Am4-see Figures 3 and S3). The free-energy barrier preventing complete translocation of NH 4 þ is astronomical (;100 À 150 RT, see Figure S3), and is located near the center of the pore around the central NH 3 binding site, Am3.

Free-Energetic Analysis of NH 3 Translocation
We also calculated the PMF profile associated with translocation of a single NH 3 molecule from the hydrated periplasmic vestibule to the cytoplasmic vestibule (see Figure  4). The profile produces three minima (sites) that were found in exactly the NH 3 positions of Am2 (z ¼ 0.17 nm), Am3 (z ¼ À0.27 nm), and Am4 (z ¼ À0.68 nm) from diffraction studies [17], thus we label them accordingly. Small ''oscillations'' in the profile are observed at z ' 0.4 À 1.0 nm corresponding to NH 3 passage across the phenyl groups, F107 and F215, forming the ''ceiling'' of the pore. A similar feature is observed on the cytoplasmic end of the pore just after site Am4 at z ' À1.1 nm corresponding to NH 3 passage across the phenyl side chain of F31, forming the ''floor'' of the hydrophobic pore (lumen). Residues H168 and H318 are seen to form stabilizing interactions with NH 3 at each binding site (Figure 4, bottom). The most stable site at single NH 3 occupancy is seen to be Am4, at À3.46 RT units. The barrier for bringing NH 3 from Am4 to one of the NH 4 þ sites (Am1 or Am5) is drastically lower (;9 À 16 RT, see Figure 4) than the barrier for bringing NH 4 þ to any one of the NH 3 sites (Am2, Am3, or Am4, see Figure S3). The structure of the pore was seen to be structurally invariant regardless of whether a single NH 3 molecule occupied Am2, Am3, or Am4 ( Figure 5). Unlike other channels such as the K þ channel, whose conformation and activity is highly dependent upon the occupancy of permeant species [30,31], the AmtB channel architecture is stiff enough to maintain the structure of all three binding sites (Am1, Am2, and Am3), regardless of the occupancy of the pore. This makes our analysis of the pore convenient, because it means that the positions of the binding sites do not change even at higher occupancies. Thus, we may rely on a 1-D PMF ( Figure  4) to enumerate the possible NH 3 loading states.
Apparent pK a of NH 4 þ as a Function of the Transport Axis Since we know the free energy as a function of the transport axis in the case of either the protonated or unprotonated form of Am, it is possible to glean the apparent pK a of NH 4 þ (pK a (z)) at all points along the transport pathway.
This will be particularly useful in identifying precisely where its deprotonation occurs, and may even tell us how. The freeenergy cycle in Figure 6A implies that we may calculate pK a (z) using Equation 1 (assuming RT units), pK a ðzÞ ¼ pK a ðbulkÞ À ðDG bulk dp À DG z dp Þ 1:2303 where pK a (bulk) is the pK a of NH 4 þ in the bulk (known experimentally to have a value of 9.25), far away from the membrane surface, DG bulk dp , is the free energy for the alchemical reaction NH 4 þ ! NH 3 in the bulk, far from the membrane, and DG z dp is the free energy for the same reaction at a particular position, z, along the transport axis (''dp'' stands for ''deprotonation''). We used the method of thermodynamic integration to obtain the necessary freeenergy values (shown in Table 1) leading to the appropriate difference curve from the PMF profiles in Figures 3 and 4 (see Materials and Methods) yielding pK a (z).
An explanation of the meaning of ''apparent'' pK a (or Àlog K a ) may be useful to the unacquainted reader. For a position, z, along the transport axis, the equilibrium constant for the proton dissociation event, where the square brackets denote populations (or concentrations). This quantity is not particularly useful experimentally (due to the ill-defined nature of the local pH or [H þ ] z ) or in conventional MD simulations (due to the inherent design of water models, which represent only  Figure 2. Am (de)protonation regions are marked with green vertical bars, whose widths give the uncertainty in the regions, and NH 4 þ binding sites are marked with gray vertical bars. The PMF profile for NH 3 (solid line) has been shifted such that its value coincides with that of NH 4 þ (dashed black line) at the equivalence ((de)protonation) regions (at the center of each green bar). Thus, the free-energy barrier for bringing NH 4 þ from site Am2 to the periplasmic equivalence point is ;15 RT units. NH 3 binding sites are marked Am2 (z ¼ 0.17 nm), Am3 (z ¼ À0.27 nm), and Am4 (z ¼ À0.68 nm). The free-energy barrier for bringing NH 3 (assuming single occupancy) from site Am4 to the cytoplasmic equivalence point is ;10 RT units. The bound state of NH 3 at each of the three sites is shown in molecular detail (bottom panel). All specific interactions within each site are described in the text. doi:10.1371/journal.pcbi.0030022.g004 neutral pH). Thus, we turn to the ''apparent'' equilibrium constant, þ ] z to describe the favorability for Am to exist in its protonated or unprotonated form. This allows us to use the fixed value for pH for the bulk, which is well-represented by our force field, pH bulk ¼ 7. With this definition, a (de)protonation event will occur at the socalled ''equivalence point,'' when the apparent pK a ¼ pH bulk , or, in other words, when pK a ¼ 7. One must note that although, normally, the auto-ionization of water provides a range of possible values for the pK a of 0-14 in aqueous solution, in a nonaqueous environment (say, in the absence of water or in an environment where water is nonbulk-like) this range does not apply.
The resulting NH 4 þ pK a profile is shown in Figure 6B.
Interestingly, the apparent pK a is shifted upward by ;5 units near the membrane surface and around sites Am1 and Am5. More precisely, using values from Table 1, the shift in pK a is þ4(62) units at site Am1 (in agreement with a previous calculation [22]) and þ5(63) units at site Am5. This upward shift, disfavoring deprotonation, is not unexpected for several reasons. First, any region that stabilizes or binds NH 4 þ will necessarily foster the protonated form of Am. This is true, not only for the binding sites, Am1 and Am5, but also for the membrane surface itself (see Figures 2 and S2), which is seen to bind cations strongly-a result commonly seen at the membrane/electrolyte interface [26][27][28][29]. In addition, water has been shown to be highly polarized near the membrane surface and shows very different electrostatic properties than seen in the bulk [32]. The above arguments rationalize the upward shift in pK a near the membrane surface and stabilization of Am in its protonated form at the Am1 and Am5 binding sites. Note that at some regions in the curve the apparent pK a is seen to have values larger than 14, which violates limits implied by the auto-ionization of water. At these regions, DG z dp .135 RT units, and deprotonation may be considered impossible. In addition, it is worth noting that the uncertainty accumulated in the thermodynamic integra-tion calculations (Table 1) and in shifting the PMF profiles leads to an uncertainty in the pK a shift of ;3 units.
In the center of the hydrophobic pore containing sites Am2, Am3, and Am4 ( Figure 6B, vertical red bars), we see that the pK a of NH 4 þ becomes negative, again violating limits implied by water auto-ionization. In this region, DG z dp , 103 RT units, which implies that it is impossible for Am to exist in its protonated form. In fact, the pK a profile indicates that Am must exist as NH 3 well before reaching site Am2 on the periplasmic end, and well after leaving site Am4 on the cytoplasmic end of the pore (see Figure 6B, in between the vertical green bars). This result contrasts with qualitative analyses of static X-ray [17] or dynamic [23] channel structures, which posit that NH 4 þ can reach site Am2 before becoming deprotonated.

The Periplasmic NH 4 þ Deprotonation Site
The equivalence point is where one would expect to find Am in the form of NH 3 and NH 4 þ with equal probability, and is marked by the two green vertical bars in Figure 6B (we also show the bars in Figures 3 and 4). At neutral bulk pH, this is where deprotonation of NH 4 þ occurs. The width of each bar covers the range of certainty (;3 pH units) for which the apparent pK a is equal to the bulk pH, or, in our case, where pK a (z) ¼ 7. Since the pK a (z) profile is extremely steep in this region, we may narrow down the deprotonation regions to very specific locations for both the periplasmic and cytoplasmic sides, despite the uncertainty in the pK a profile.
On the periplasmic end of the channel, the site of deprotonation is immediately beneath site Am1 at z ¼ 0.77 6 0.08 nm (with respect to the center of mass of the trimer backbone). This site coincides exactly with the position of the phenyl side chain of F107 (see Figure 6C). As NH 4 þ travels from site Am1 toward the cytoplasm, its pK a drops so steeply that it must become NH 3 before passing the second phenyl group, F215. Based on Figure 4, the process of bringing NH 4 þ from site Am1 to the deprotonation region corresponds to crossing a free high-energy barrier of ;15 RT units, which supports the experimental observation of very slow Am translocation rates of 10 5 -10 8 ns/molecule [18]. If we observe the average hydration number of Am as a function of the transport axis ( Figure 6D), we can rationalize the extremely unfavorable (15 RT) barrier to translocation and drastic drop in apparent pK a . In the bulk, NH 4 þ is hydrated by six water molecules, and at site Am1 (vertical gray bar at z ; 1.07 nm) it loses two to three water molecules from its hydration shell to form hydrogen bonds with A162 and S219 ( Figure 3), but as the external electrochemical gradient pushes it toward the deprotonation site, its hydrogen bonds are stripped away such that eventually the only available acceptors are one to two water molecules, and the carbonyl oxygens of A162 and F215. The transition state for the NH 4 þ to NH 3 transformation at the deprotonation site may be most accurately represented (given our classical molecular description) by Figure 7, where the F107 phenyl group is seen to allow passage by rotating in response to the presence of the cation. In this configuration, NH 4 þ is stripped to three hydrogen bonds with one water molecule and the carbonyl groups of A162 and F215, and must release one of its protons in order to favorably continue toward the cytoplasm as NH 3 . Previous works have suggested that the proton is donated to A162 [23] and, eventually, to  . Apparent pK a and Hydration of Am Along the Transport Axis (A) Free-energy cycle used to determine pK a of NH 4 þ as a function of the transport axis. DG bulk dp is the free energy for the alchemical reaction NH 4 þ ! NH 3 in the bulk, far from the membrane, and DG z dp is the free energy for the same reaction at a particular position, z, along the transport axis. DG z NH4 is the free energy to bring NH 4 þ from the bulk to a particular position, z, and DG z NH3 is the analogous quantity for NH 3 . þ binding sites, red bars represent NH 3 binding sites, and green bars represent (de)protonation regions. The periplasmic (de)protonation region is at z ¼ 0.77 6 0.08 nm, and the cytoplasmic (de)protonation region is at z ¼ À1.13 6 0.07 nm. See the text for an explanation of the uncertainties (width of the green bars) in these regions. (C) Synopsis of binding sites and deprotonation regions generated by multiple structural alignment of AmtB channel structures singly bound at each binding site (a total of five structures taken from umbrella sampling simulations) and the X-ray structure [17]. Some transmembrane helices have been removed from the protein (orange) to allow visibility of the binding sites. The positions of the X-ray sites are shown as red spheres, and Am binding sites derived from our analysis are shown as either blue spheres (in the case of NH 3 ) or blue spheres with white hydrogen atoms (in the case of NH 4 þ ). Any observed differences in the positions of the sites are within fluctuations. (De)protonation regions are marked in green. The periplasmic and cytoplasmic (de)protonation regions coincide with the phenyl groups of F107 and F31, respectively. (D) Hydration of NH 3 (diamonds with dashed line) and NH 4 þ (circles with solid line) as a function of the transport axis derived by combining statistics from free MD production runs and umbrella sampling runs. The hydration number is defined as the average number of water molecules falling within the first coordination shell of Am. The first coordination shell of Am was taken to be the position of the first minimum in the N Am À O water pair correlation function (R coord ¼ 0.362 nm, see Figure S4). doi:10.1371/journal.pcbi.0030022.g006 D160 [20,23]. However, at this deprotonation site, NH 4 þ has access to periplasmic vestibular water, which connects to the bulk in a continuous fashion ( Figures 6D and 7). Thus, it is much more likely to pass its proton directly to water, so that the proton may escape directly to the periplasm in the form of hydronium (cutting out the ''middle man''). This conclusion is strongly supported by (1) the fact that water ionizes more easily than a carbonyl group; (2) the fact that at the transition state ( Figure 7) the carboxylate of D160 is persistently engaged in hydrogen bonds with the hydroxyl and amide nitrogen of T165, and the amide nitrogens of G163 and G164; and (3) the finding that these strong interactions shift the pK a of the D160 sidechain downward by ;9 pH units, making its protonation effectively impossible [22]. Therefore, it appears that the role of D160 is not to accept a proton, but, as mentioned in structural analyses [17,18], to force recruitment of cations from the bulk by orienting periplasmic M4-M5 loop carbonyl groups, and by electrostatically stabilizing NH 4 þ at site Am1 [22]. We again emphasize that these roles for D160 are likely not mutually exclusive.
It is also interesting to note that Figure 6D shows that the average hydration number of Am at the periplasmic deprotonation region (vertical green bar at 0.77 nm) drops immediately from 1-2 to 0 when NH 4 þ becomes NH 3 , as if this transition immediately allows Am to shed its final water to favorably enter the hydrophobic pore. In this sense, even though the reaction occurs in a condensed phase, it is much like the reaction that would occur when transferring Am from aqueous solution to a gas phase. The process of bringing Am from site Am1 (as NH 4 þ ) to site Am2 (as NH 3 ) is summarized in Figure 8A.

The Cytoplasmic NH 3 Reprotonation Site
After NH 4 þ is dehydrated and deprotonated at the periplasmic deprotonation region (near F107), it continues as NH 3 , completely dehydrated ( Figure 6D), and moves across sites Am2, Am3, and Am4. The consistent slope in the PMF (Figure 4), although all three NH 3 binding sites provide stable sites, pushes NH 3 toward the most stable site, Am4, on the periplasmic end of the hydrophobic lumen. Here, further passage is hindered by the phenyl group of F31. This sidechain, forming the ''floor'' of the hydrophobic lumen provides an ;10 RT free-energy barrier (see Figure 4) for NH 3 reaching the cytoplasmic (de)protonation region at z ¼ À1.13 6 0.07 nm (see the corresponding green bar in Figure  6B-6D). This region nearly coincides with the position of the phenyl group of F31 ( Figure 6C). The nature of the barrier to reprotonate NH 3 is seen to be similar in origin to that for deprotonation of NH 4 þ near F107.
Namely, it requires work to force NH 3 toward the wateraccessible cytoplasmic vestibule just beneath the F31 phenyl ''floor'' of the lumen. Presumably, larger NH 3 luminal occupancy (i.e., more than 1 NH 3 ) would lower the barrier to reprotonation due to mutual destabilization of NH 3 molecules (we discuss this later). The protonation of NH 3 to form NH 4 þ at the reprotonation region (green bar in Figure   6D) is coupled to a gain in hydration from 0 (coordinating NH 3 ) to 2-3 (coordinating NH 4 þ ) water molecules and the gain of one hydrogen bond to the hydroxyl oxygen of S263 (see Figure 6D and Figure 8B3). Thus, again, the only plausible species involved in the protonation event is water. The process of bringing Am from site Am4 (as NH 3 ) to site Am5 (as NH 4 þ ) is summarized in Figure 8B.

Selective Permeability: The Role of Phenylalanine and the (Dehydrated) Hydrophobic Lumen
Given the results presented thus far, it is clear that phenyl side chains play a key role at both the periplasmic and cytoplasmic (de)protonation regions. At the periplasmic (de)protonation region, the presence of the phenyl sidechain All values from NH 4 þ ! NH 3 reactions were subsequently used in pK a profile determination (see Materials and Methods). The NH 3 ! 0 reactions represent the alchemical transformation of an NH 3 molecule to a ''null'' molecule restrained to its site (represented by ''0'') which has no interaction with the rest of the system. These reactions were used to complete the map in Figure 9A (see Materials and Methods and Figure S5). b The AmtB channel binding site at which the Am species was located during the alchemical transformation.  from F107, though it rotates to allow passage, provides a constriction in the Am pathway such that most of its hydrogen bond acceptors (including water) are removed ( Figure 7). Given the sequence of events observed in Figure  8A, it appears that once NH 4 þ donates a proton to neighboring water and moves slightly farther into the lumen, the F107 side chain guards against hydration while the second phenyl side chain from F215 rotates to allow passage of NH 3 . Together, F107 and F215 act as a ''double doorway,'' aiding to ensure dehydration of the passing Am species. In a similar manner, at the floor of the lumen, F31 guards against hydration from the cytoplasm ( Figure 8B). During 27.5 ns of equilibration and 25 ns of production MD with no intraluminal Am species, neither water nor NH 4 þ was seen to spontaneously enter the hydrophobic lumen (see Figure S4). During umbrella sampling simulations where NH 4 þ was forced to sample the lumen, it was seen to carry along anywhere from zero to three water molecules (compared with six in the bulk, see Figure 6D). However, as mentioned before, the free-energy barrier for NH 4 þ entering the lumen from either the periplasm or cytoplasm is prohibitive (.15 RT to reach the deprotonation region and .50 RT to reach the NH 3 sites Am2, Am3, or Am4, see Figures 3 and S3). On the other hand, NH 3 is quite ''happy'' to be completely dehydrated as it occupies the lumen (see Figures 4 and 6D). Therefore, we can conclude that the NH 4 þ ion's ability to deprotonate to form NH 3 imparts Am with the ability to permeate the AmtB channel in a selective fashion.  That is, permeation of Am or methylated Am will be greatly preferred over other aqueous inorganic ionic species because deprotonated Am (or methylated Am) can favorably occupy the hydrophobic lumen. Based on our calculations, we can expect prohibitive free energy barriers (.50 RT) for permanently charged species to enter the lumen. Thus, without a doubt, this is the most robust rationale for the AmtB channel's display of selective permeability for Am over all permanently charged species (e.g., Family IA and IIA cations such as K þ or Na þ ), and is the most plausible source of the ''Am-sensing'' property of this channel [24]. With that being the case, the Phe residues, F107, F215, and F31, play a role in selective permeation, not by providing a larger relative stability at the periplasmic/cytoplasmic binding sites (Am1/Am5) for NH 4 þ over other cations, but by aiding in the removal of water at the deprotonation region of the transport pathway. Whether other ions can competitively inhibit binding, and thus recruitment, of NH 4 þ to the periplasmic site, Am1, is a different question altogether (aside from selective permeability). However, this appears not to be an issue since experimental studies do not observe significant competitive inhibition of methyl-Am uptake by alkali cations in Amt/MEP/ Rh transporters [6,33]. Recent calculations [22] reported slight selectivity of site Am1 for NH 4 þ over monovalent Family IA cations, rationalizing experimental observations of mild Am uptake inhibition. With only mild selectivity at site Am1, it appears that the hydrophobic lumen is the most Amselective feature of the channel. The finding that Am sheds its water before entering the hydrophobic lumen as NH 3 , where it remains dehydrated, contrasts with a recent computational study by Nygaard et al., which reported that the interior lumen of the hydrophobic pore is hydrated [23]. In that study, two 10-ns simulations, each with different neutral tautomeric states of the intraluminal His pair, H168 and H316, both demonstrated hydration of the lumen all the way up to the periplasmic phenyl ''stack'' provided by F107 and F215. The water entered the lumen in the form of a highly ordered, hydrogen bonded ''water wire,'' a very commonly observed arrangement in hydrophobic pores [34,35], which entered the pore via the cytoplasmic vestibule [23]. This observation led the authors to suggest that proton conduction by single-file water might provide a mechanism for H þ exchange with H168 and H316.
In contrast to the study of Nygaard, our MD simulations never displayed this hydration event in any one of the three channels in the AmtB trimer, despite 52.5 ns of dynamics (combining both equilibration and production runs). It is possible that a difference in choice of force field (see Materials and Methods) may have been the origin of these contrasting results. However, since the study of Nygaard et al. included only Cl À counterions so as to achieve a net neutral system charge (no electrolyte) [23], we find it more probable that the main cause for our different observations is inclusion of NH 4 Cl electrolyte in our simulated system (measured to be 158 6 5 mM in the bulk, see Figure 2B). We say this because past work has shown that simulated systems utilizing 3-D periodicity with full Ewald sums for treatment of electrostatic interactions and comprising membranebound proteins possessing a significant dipole moment (greater than ;300 Debye) can display unexpected ''water-ordering'' artifacts [36,37]. This water-ordering phenomenon is constant in the bulk, can propagate all the way into the pores of membrane-bound channels, and is caused by the infinite array of persistent dipoles (caused by the membranebound protein in the central simulation cell) represented by a fully periodic system in 3-D. Normally, the largest dipole component of a membrane-bound protein will fall along the membrane normal (or transport axis). Thus, one way to remove the effect due to the electric field produced by the infinite, 3-D periodic dipolar array is to remove the periodicity in this direction (the membrane normal is usually taken to be the z-dimension), and to simulate in a slab (effectively 2-D or xy-) geometry [36]. There are other ways to avoid this artifact. For example, one can add electrolyte to the system (as we have done) so that it bears the dielectric response to the protein, rather than the water [37]. One may also choose to use a cutoff treatment of electrostatic interactions [37], although other artifacts are associated with this method. The addition of counterions alone, though it is important to assure a neutral system charge, does not appear to be enough to avoid the effect [36,37]. Nor does the addition of generous layers of water in bulk in an attempt to screen the effect [37].
The dipole moment of the AmtB trimer is several times larger (;2,000 Debye, see Figure S1A) than the membrane proteins for which this phenomenon was first observed, and thus is an excellent candidate for providing an electric field to which the surrounding simulated system must respond. In our system, we observe that the dipole moment of the electrolyte precisely counteracts the dipole of the protein (Figure S1A), leaving a net zero dipole moment for the entire system ( Figure S1B). Without the electrolyte, we can expect the water to bear the majority of the burden of the dielectric response to the protein.
The next obvious question is: Given that the net force of an electric field on a dipolar molecule is zero, why might this cause water to enter the hydrophobic lumen of AmtB? Recently, much attention has been given to the occupancy of water in hydrophobic spaces provided by macromolecular assemblies, including narrow hydrophobic pores, cavities, and the hydrophobic membrane interior [34,35,38,39]. It has been shown that, indeed, electric fields directed parallel to the axis of a narrow hydrophobic pore can make a water-filled state favorable, whereas in the absence of a field, it would be unfavorable [35]. In addition, electric fields directed along the bilayer normal can cause water to penetrate the hydrophobic portion of the lipid membrane itself and form water-filled pores-a process known as electroporation [38]. In either case, the electric field induces a preferred water orientation, which increases the probability of observing ordered water defects either at the membrane surface (in the case of membrane electroporation) or at the hydrophobic pore entrance (in the case of pore hydration). This appears to increase the probability that water enters the pore. It is also known that hydrophobic pore hydration occurs in an ''all-ornone'' fashion, well-represented by two-state models [34,35]. Once water enters a hydrophobic pore, we can expect that it will become filled. That is, one would assume, unless there is some feature that blocks full permeation. It would appear that the phenyl groups of F107 and F215 provide such a ''blocking'' feature.

Free Energy and Dynamics of Multiply NH 3 -Occupied Luminal States
The three NH 3 binding sites within the hydrophobic lumen proved to be highly disordered in the X-ray structure [17], leading to uncertainty in understanding whether the sites were occupied alternately with one another during the conductance mechanism, or partially occupied such that at high-bulk Am concentrations they might be fully occupied. According to the PMF for NH 3 at single occupancy (Figure 4), it appears that immediately after the deprotonation event, NH 3 is most likely to jump across sites Am2 and Am3 to site Am4. There it must either overcome a barrier (;37 kJ/mol or ;15 RT) to reprotonate and escape back to the periplasm, or overcome a barrier (;25 kJ/mol or ;10 RT) to reprotonate and escape to the cytoplasm. Multiple occupancy of the lumen will undoubtedly lower these barriers.
We utilized additional thermodynamic integration calculations (see Table 1 and Materials and Methods) to draw a state map for the possible NH 3 -occupied states of the lumen, shown in Figure 9A. For brevity, we adopt a vector representation for the occupied state of the three-site lumen, (O 2 ,O 3 ,O 4 ), where O 2 , O 3 , and O 4 are binary digits indicating whether Am2, Am3, or Am4, respectively, are occupied (O i ¼ 1) or unoccupied (O i ¼ 0). Since NH 3 most favorably occupies site Am4 at single occupancy, (i.e., state 001), the map suggests that when a second NH 3 enters, it will most likely occupy site Am2 (i.e., state 101). Figure 9A suggests that NH 3 binding at site Am4 when site Am2 is occupied (the reaction 100 ! 101) is 7 kJ/mol less-favorable than binding at site Am4 when no other site is occupied (the reaction 000 ! 001). If we assume that the barrier for NH 3 escape from Am4 to the cytoplasmic protonation region remains the same (we will refer to this as the ''barrier assumption'') as calculated for single occupancy (Figure 4), then it would cost less free energy (;18 kJ/mol or 7RT) for NH 3 to move on to the cytoplasm from state 101. If NH 3 were to dynamically move to a state where Am3 and Am4 are simultaneously occupied (state 011), we find that NH 3 binding at site Am4 when site Am3 is occupied (010 ! 011) is 28 kJ/mol less-favorable than binding at site Am4 when no other site is occupied (000 ! 001). And with the above barrier assumption, the translocation of NH 3 to the deprotonation region from site Am4 would become favorable, requiring ;À3 kJ/mol or ÀRT. Needless to say, at full occupancy (state 111), the same approximation suggests that passage of NH 3 from site Am4 to the deprotonation region would be even more favorable (;À10 kJ/mol or À4RT).
Thus, our results suggest that doubly-occupied NH 3 dynamics in the lumen can force the escape of NH 3 at site Am4 to the cytoplasmic deprotonation region for successive binding to site Am5 and escape to the cytoplasm. The state most likely preceding the escape of NH 3 to the deprotonation region is one where the lumen is occupied at sites Am3 and Am4 (state 011). We also note, however, that the free energy to attain full occupancy (the reaction 011 ! 111) is favorable ;À9 kJ/mol or À4 RT. In fact, this reaction is more favorable than the escape of NH 3 to the cytoplasm. Therefore, it is also possible for the fully occupied state, 111, to precede escape of NH 3 to the deprotonation region.
Though, above, we have only discussed the process of NH 3 translocation toward the cytoplasm, our results do not preclude translocation toward the periplasm, which is wellknown to occur [40]. In fact, at full occupancy (state 111), it becomes more favorable (by ;À13 kJ/mol or À5 RT) for NH 3 to escape from site Am2 to the periplasmic deprotonation region (i.e., the reaction 111 ! 011) than the same reaction at zero occupancy (i.e., the reaction 100 ! 000). Thus, our results are also consistent with the facilitation of gradientdependent bidirectional NH 3 transport.
In Figure 9B, we show a brief demonstration of NH 3 dynamics for the fully occupied state of the lumen. Although the channel is fully occupied, it appears that the qualitative trend in binding-site favorability represented by the PMF of Figure 4 for single occupancy persists. The interactions of NH 3 with the lumen display a strong tendency to force NH 3 toward the cytoplasmic binding site, Am4. At times, the position of Am4 can be transiently occupied by two NH 3 molecules (see Figure 9B and 9C). We can expect slight differences in the free energetics and dynamics of NH 3 in the lumen depending on the tautomeric form of the intraluminal histidine residues [23]. The tautomeric configuration repre-  Figure 4, red values were determined algebraically by completing the embedded thermodynamic cycles, and blue values were determined using thermodynamic integration techniques (see Figure S5 and surrounding discussion). (B) A demonstration of NH 3 dynamics (1 ns) in the lumen. Horizontal red bars represent sites (listed in decreasing value of z) Am2, Am3, and Am4, respectively. For the first 300 ps of the simulation, the NH 3 molecules are restrained (in the z-dimension only) to their binding positions. When the restraint is released, NH 3 is seen to behave dynamically in the lumen, and there is a strong tendency for NH 3 to move toward the site closest to the cytoplasmic end of the channel, site Am4. (C) A snapshot illustrating the tendency for NH 3 to migrate toward site Am4. Here three NH 3 molecules have migrated from a state where Am2, Am3, and Am4 are occupied, to a state where one molecule occupies Am3 and two molecules ''fight'' to occupy site Am4. doi:10.1371/journal.pcbi.0030022.g009 sented by our simulations appears to favor translocation toward the cytoplasm, while the study of Nygaard et al. showed that the alternate form for H168 and H318 does not appear to favor translocation in either direction. It is possible that, together, these different forms may be used by the channel to facilitate translocation toward the cytoplasm or toward the periplasm, depending on the electrochemical gradient conditions. Thus, our results suggest that a periplasmic electrochemical gradient will push NH 4 þ from site Am1 to be deprotonated at the deprotonation region; the hydrophobic NH 3 molecules begin to fill up the lumen, resulting in mutual destabilization or a buildup of ''pressure,'' which in turn allows NH 3 reprotonation to become favorable. This means that the rate-limiting step of Am permeation is the NH 4 þ deprotonation step at the periplasmic end of the channel, representing an ;15 RT freeenergy barrier (see Figures 3 and 4) that must occur at the position of the F107 phenyl ring, before proceeding into the lumen, just underneath the F215 phenyl ring.

Conclusions
Our results show that the global electrostatic nature of the AmtB channel inhibits anion binding to the periplasmic membrane and plays a large role in the recruitment of cations to the outer vestibule. PMF analysis verifies all crystallographically observed [17]  passes and is stripped of most of its hydrogen bonds. At this point, only the backbone carbonyl groups of A162 and F215 and a water molecule serve as hydrogen bond acceptors. Given (1) that water ionizes more easily than a carbonyl group, (2) that the carboxylate of D160 is persistently engaged in hydrogen bonds with residues of the M4-M5 loop, and (3) that the apparent pK a of D160 is shifted downwardly by ;9 pH units due to its stabilized interactions [22], water is the only plausible candidate for accepting a proton from NH 4 þ at the deprotonation site. The proton has full access to the periplasmic bulk, and must therefore escape in the form of hydronium. Thus, the Am deprotonation mechanism of AmtB is due to stripping NH 4 þ hydration to a critical point, much like the deprotonation event that one would observe in transferring Am from a hydrated state to a gas state. Recently, a mutagenesis study showed that the twin luminal His residues (H168 and H318) are essential for substrate conductance [41]. Javelle et al. offered two explanations for this result: (1) these His residues serve as proton acceptors for entering Am such that they may move across the lumen as NH 3 , or (2) efficient conductance might require a ''narrow mainly hydrophobic pore with a few precisely oriented hydrogen bond acceptor or donor functions for weak, stabilizing interactions with water/ammonia that still permit rapid diffusion.'' Our results suggest that the latter explanation (2) applies.
The need of NH 4 þ for hydration, or some other form of hydrogen-bond stabilization, also manifests itself in our PMF analysis as an insurmountable free-energy barrier (.50 RT), preventing passage across the phenyl group of F107 and approaching any site within the lumen. This result contrasts with existing hypotheses [17,23] that NH 4 þ can occupy the intraluminal NH 3 sites, Am2, or Am3. Since Am has the rare ability to deprotonate, it may enter the lumen as NH 3 . This is the origin of the Am-sensing ability of AmtB. Since most cations, such as Family IA and IIA ions, are permanently charged, their passage is expected to be disfavored, with a free-energy barrier of more than 50 RT.
Due to the nature of the deprotonation mechanism we describe, fully ordered hydration of the lumen would inhibit the Am-sensing ability of AmtB, providing the possibility for luminal Am to exist in its protonated form. Although, in contrast to other simulations of AmtB [23], the MD simulations of our work do not show any ordered luminal hydration (i.e., in the form of a water wire), they do not rule out the possibility of a few water molecules ''sneaking'' into the lumen. Such events would not be expected to inhibit NH 4 þ deprotonation (and, thus, AmtB's sensing capability).
Recent diffraction studies [17,41] clearly show electron density at the three luminal NH 3 binding sites (as shown in Figure 4), but are unable to determine whether these sites are occupied by NH 3 or water molecules. The diffraction study of Khademi et al. [17] could not find density corresponding to these sites for AmtB structures crystallized in the absence of 25mM AmSO 4 . This suggests that these peaks are mostly due to NH 3 rather than water. Upon the introduction of charged species to the hydrophobic lumen, we might expect water to enter. For example, in the very unlikely event that NH 4 þ enters the lumen, we show that it can carry water with it (see Figure 6D). Recently, an X-ray diffraction study of an H168D mutant of AmtB showed electron density inside the lumen that might be attributable to water [41]. It was deemed most likely that this was due to the negative charge of the Asp side chain protruding into the lumen.
Given that hydration of narrow hydrophobic pores in the absence of an electric field is known to be unfavorable [35], the absence of ordered luminal hydration is most likely attributable to the presence of NH 4 Cl electrolyte in our simulated system, which, by responding to the macro-dipole of AmtB in the membrane, neutralizes the resulting electric field acting parallel to the transport axis. Similar ''neutralizing'' electrolytic responses to dipolar membrane proteins have also been observed in other work [37]. It appears that the periplasmic phenyl groups of F107 and F215 and the cytoplasmic phenyl group of F31 are major contributors to the maintenance of the dehydrated state of the AmtB lumen. This is most likely the reason why our NH 4 þ pK a calculations highlight the positions of F107 and F31 as deprotonation ''landmarks.'' Once Am occupies the lumen as NH 3 , in the singly occupied state, Am4 near the cytoplasm is the most stable site. As more hydrophobic NH 3 is added to the lumen by means of NH 4 þ deprotonation from the cytoplasm, it may favorably become reprotonated from either a doubly occupied state (where Am3 and Am4 are simultaneously occupied) or a triply occupied state. Favorability for reaching the cytoplasmic reprotonation site is thus attained by a buildup of NH 3 pressure in the lumen. In addition, the dynamics of NH 3 within the lumen supports flow toward the floor of the lumen at the cytoplasmic site, Am4. Since, at the floor of the lumen, reprotonation of NH 3 occurs favorably at multiply occupied states, the rate-limiting step for Am permeation must be the deprotonation step at the periplasm, where NH 4 þ becomes dehydrated. The free-energy barrier for this step is ;15 RT based on our PMF calculations. Though our study addresses mostly the translocation of Am toward the cytoplasm, the results we present are not inconsistent with the expected capability of bidirectional Am translocation. A previous study has found slight differences in the dynamics of NH 3 in the lumen depending on the tautomeric form of the intraluminal histidine residues [23]. It is possible that, together, these different forms may be used by the channel to facilitate translocation toward the cytoplasm or toward the periplasm, depending on the electrochemical gradient.
Based on our apparent pK a calculations, reprotonation of NH 3 at the cytoplasmic end of the lumen was shown to occur upon hydration immediately after passing the phenyl group of F31. No plausible donor other than water exists at the point where NH 3 reaches the reprotonation site. Thus, reprotonation at the cytoplasmic end of the channel must occur in the reverse manner to deprotonation at the periplasmic end. Upon reprotonation, Am binds as NH 4 þ to site Am5 before escaping to the cytoplasm. Thus, the role of the newly discovered site, Am5, is apparently to reduce the barrier to reprotonation by promoting hydration directly beneath F31 and by stabilizing NH 4 þ such that it does not escape back into the lumen as NH 3 .

Materials and Methods
Setup, equilibration, and production simulations of the AmtB system. All MD simulations were performed using the GROMACS package [42,43], with the OPLS force field for protein, ions, and NH 3 [44,45], the OPLS lipid parameters developed by Smondyrev and Berkowitz [46], and the SPC water model [47]. Our choice in force field parameters was guided by the fact that the OPLS force field, when used with the SPC water model, has been shown to yield equivalent (or better) results than the standard TIP3P water model when determining solvation free energies of amino acid side-chain analogs [48].
Production runs used a time step of 2.0 fs. Configurations of the system were saved every 1 ps for later analysis. Periodic boundary conditions were applied in all three dimensions. The LINCS algorithm was used to constrain all bonds in the systems [49]. Longrange electrostatics were handled using the PME algorithm [50], with a real-space cutoff of 0.9 nm. Other nonbonded interactions were truncated at 1.2 nm. The temperature was maintained at 298 K using the Nose-Hoover scheme with an oscillatory relaxation period of 2.0 ps. The pressure was maintained at 1.0 atm using the Parrinello-Rahman coupling scheme [51,52] with a barostat time constant of 2.0 ps. The rectangular simulation box was allowed to scale in size semiisotropically to maintain pressure. Analysis of trajectories was performed using a combination of GROMACS analysis utilities and locally written scripts and programs.
The initial configuration of the AmtB trimer system was generated by first equilibrating a hydrated POPC bilayer containing 450 lipids (225 lipids per leaflet) and 29,889 water molecules for 4.0 ns. The area per lipid molecule was seen to converge to ;0.62 nm 2 . The X-ray structure of Khademi et al. (PDB 1u7g) was used as the initial structure of the AmtB trimer in our simulations [17]. All residues were mutated such that a wild-type representation of the protein was obtained. In addition, the missing residues, A1 and P2, were added to each monomer to yield the full protein. All crystallographic watermolecule positions that coincided with the transmembrane portion of the trimer were stripped from the structure except for those deemed by Khademi et al. to be in ''special positions.'' We also removed NH 3 positions from the interior binding sites (Am2, Am3, and Am4), to be replaced later for other calculations, leaving only the NH 4 þ position at site Am1 for each monomer. Hydrogen atoms were then added to the protein, NH 4 þ molecules, and crystallographic water molecules. All residues were assigned their standard (neutral pH) protonation states, and the tautomeric forms of each His residue were assigned using (geometric) hydrogen-bonding criteria except for the intraluminal H168 and H316 residues, which were shown to share a proton, and were thus assigned the tautomeric form suggested by Khademi et al. [17] and shown in Figure 4. The resulting structure used for subsequent equilibration contained three AmtB monomers, three NH 4 þ ions, and 486 water molecules. A previous computational study [23] supported the tautomeric forms of H168 and H316 used in our simulations (shown in Figure 4) for the case where water is absent from the lumen. Since, for all of our simulations, water was not seen to enter the lumen, it would appear that our choice of tautomeric state is reasonable. It should be noted, however, that any conformational event observed involving these His residues during simulation (including umbrella sampling simulations) can only represent one of a few possibilities, since their protonation states must be permanently assigned. Of course, this is true of any simulation employing conventional force fields.
The trimeric structure described above was inserted into the equilibrated POPC bilayer, and all overlapping lipid and water molecules in the hydrated bilayer were removed. Water molecules in the bulk were randomly replaced by 117 Cl À and 108 NH 4 þ . The net charge on each AmtB monomer was 2e, thus with 117 Cl À ions, and a total of 111 NH 4 þ ions (three coming from the X-ray structure), the net charge of the system was neutralized. The resulting initial system contained the AmtB trimer, 117 Cl À ions, 111 NH 4 þ ions, 291 POPC molecules, and 28,099 water molecules (a total of 117,276 atomic sites).
The system was minimized using the method of steepest descent, and simulated for 1.5 ns using Berendsen temperature and pressure coupling [53] while imposing restraints on the protein to allow membrane, ions, and water to relax around the protein. Restraints on the protein were then released, and the system was equilibrated further for 1 ns, again using Berendsen coupling. We then switched to Nose-Hoover temperature coupling and Parrinello-Rahman pressure coupling, and thermalized the system further, for 27.5 ns. We show the time evolution of the center of mass of both Cl À and NH 4 þ for each leaflet of the biomembrane system ( Figure 10A). It is seen that this quantity is well-converged after the first 15 ns of the equilibration run. The protein root mean square deviation (RMSD) was seen to converge after the first 15 ns of total simulation time as well ( Figure 10B).
After the equilibration run, a 25-ns production run was performed. From this simulation, we derived the partial densities of various species across the simulation box (see Figures 2 and S2) and the hydration number of NH 4 þ across the simulation box (see Figures  6 and S4), to be used later, in conjunction with similar analyses of umbrella sampling trajectories. The average concentration of NH 4 Cl in the bulk as measured from the resulting partial densities was 158 6 5 mM (see Figure 2B).
After the 25-ns production run, we placed three NH 3 molecules into the luminal binding positions, determined from umbrella sampling simulations (Figure 4), within a single monomer of the final configuration. To demonstrate the dynamics of NH 3 in the lumen, we performed 300 ps of simulation with the NH 3 molecules harmonically restrained (in the z-dimension only-with a spring constant of 5,000 kJ/mol/nm 2 ) to their binding positions, and 700 ps of simulation with the NH 3 molecules unrestrained. The relevant data from these simulations is shown in Figure 9B.
PMF calculations. The final structure after 27.5 ns of equilibration was taken as the starting structure for umbrella sampling simulations [54]. All umbrella sampling windows utilized the same simulation conditions that were used for the production MD run. Am species (either NH 4 þ or NH 3 ) were restrained to z-positions along the transport axis with respect to the AmtB trimer backbone (Am was free to move in the xy-plane). Each window utilized a harmonic umbrella potential acting in the z-dimension only, with a spring constant of 3,000 kJ/mol/nm 2 . Each umbrella simulation consisted of 150 ps of equilibration, followed by 150 ps of production MD over which statistics of Am occupancy along the transport axis (z) were recorded. Statistics from umbrella window simulations were combined and unbiased using the weighted histogram analysis method (WHAM) [55] to yield the probability, P(z), that a given Am species occupies position z along the transport axis. The PMF in units of RT, where R is the gas constant and T is the temperature (RT ¼ 2.48 kJ/ mol at 298 K), is then given by the Boltzmann relation W(z) ¼ÀlnP(z).
Given the size of the system and the large barriers required for translocation of the two Am species (i.e., involving displacement of Phe sidechains and (de)hydration barriers), we took a conservative, incremental approach to building the PMF for single NH 4 þ and NH 3 translocation to ensure adequate relaxation of the system at each window simulation. In the first series of umbrella simulations, we sampled NH 4 þ along z from site Am1 (z ¼ 1.04 nm), outwardly, into the periplasmic solution (z ¼ 2.64 nm) using windows of Dz ¼ 0.04 nm (a total of 41 window simulations). Next, we sampled NH 4 þ along z from just beneath site Am1 (z ¼ 1.00 nm), inwardly, into the lumen of the channel (down to z ¼À0.20 nm) using windows of Dz ¼ 0.04 nm (a total of 31 window simulations). In both series of simulations described above, the initial configuration of the (i þ 1)th umbrella simulation was taken from the final configuration of the ith umbrella simulation. Thus, each simulation within a series was performed consecutively rather than in parallel (as is usually done) to ensure proper relaxation of the system along the Am pathway.
We then moved on to sample NH 3 along its transport path starting within the channel lumen using two additional series of simulations. The initial configuration for these series was generated by removing the three Am1-bound NH 4 þ ions and three Cl À ions from the bulk solution (to maintain net neutral system charge) in the initial structure of the NH 4 þ umbrella simulations described above. A single NH 3 molecule was placed at site Am3 (the central luminal site, at z ¼ 0.25 nm) within the AmtB monomer. We then sampled NH 3 along z from site Am3 (z ¼ 0.25 nm), outwardly, toward the periplasmic vestibule (up to z ¼ 1.35 nm) using windows of Dz ¼ 0.04 nm (a total of 41 window simulations). Next, we sampled NH 3 along z from just beneath site Am3 (z ¼ À0.29 nm), inwardly, toward the cytoplasmic vestibule (down to z ¼ À1.89 nm) using windows of Dz ¼ 0.04 nm (a total of 41 window simulations). Again, each simulation within a series was performed consecutively.
Finally, we sampled NH 4 þ translocation at the cytoplasmic end of the channel. The initial structure for these series of simulations was taken from the 27.5-ns equilibration run. The starting position of NH 4 þ was taken to be a position of NH 3 near the phenyl group of F31 generated from umbrella sampling of NH 3 translocation described above (at z ¼ À1.10 nm). NH 4 þ translocation was sampled outwardly, toward the periplasm (up to z ¼À0.02 nm) using windows of Dz ¼ 0.04 nm (a total of 28 window simulations). We also performed an inward sampling series of simulations, starting from z ¼À1.14 nm toward the cytoplasm (down to z ¼À2.70 to nm) using windows of Dz ¼ 0.04 nm (a total of 40 window simulations).
After combining and unbiasing all NH 4 þ translocation umbrella statistics, we matched the resulting PMF at the periplasmic and cytoplasmic ends with the PMF derived from NH 4 þ densities provided by the 25-ns MD production run to extend the profile into the bulk region ( Figure 3). The PMF was shifted such that its value in the bulk was zero. Both the NH 4 þ and NH 3 PMFs were smoothed using window averaging. Taking into account all trajectories used in the construction of the PMFs in this work, including all umbrella simulations and the production run, a total of 91.6 ns were simulated. Production run configurations were saved every 1 ps, and statistics from umbrella sampling trajectories were taken every timestep, yielding an analysis of a total of 1.69 3 10 7 configurations for constructing the PMFs. The total simulation time mentioned above (91.6 ns) does not take into account the equilibration run (27.5 ns), nor the thermodynamic integration calculations (eight calculations comprising 2.1 ns eachsee the discussion below), which, when combined, give the total simulation time of this study: 135.9 ns.
Alchemical transformations and pK a (z) determination. We performed two types of free-energy calculations in this work (results summarized in Table 1): (1) to determine the free energy for mutating NH 4 þ to NH 3 at particular locations along the transport axis, and (2) to determine the free energy for ''turning off'' all interactions between NH 3 and the rest of the system at particular binding sites within the lumen of the AmtB channel. For each individual free energy calculation, the system topology was varied, using a coupling parameter, k, from the initial state (k ¼ 0) to the final state (k ¼ 1). All transformations were performed over a set of 21 simulations, each of 100-ps length and carried out at different fixed values of k (i.e., k ¼ f0.00,0.05,0.10,. . .,0.95,1.00g). Data from the final 90 ps of each simulation was split into nine blocks, each 10 ps in length, thus providing nine sets of data (21 values each for h@H/@ki P,T , where H(k) is the hamiltonian of the system). These sets of data were integrated from 0 to 1 to obtain the free energy of the transformation, The nine resulting values for the free energy were averaged to obtain the final value, and an upper bound for the uncertainty in this value was taken to be the standard deviation of the sample. The type 1 calculations were used in conjunction with the PMF profiles of Figures 3, 4, and S3 to determine the apparent pK a as a function of the transport axis ( Figure 6B). In these calculations, the NH 4 þ molecule (k ¼ 0) of interest was transformed to an NH 3 molecule bonded to a dummy atom (k ¼ 1), where the dummy atom had no interactions with the remainder of the system. The pK a calculation required (see Equation 1) the free-energy difference, DG bulk dp , for the NH 4 þ ! NH 3 transformation at some position in the bulk. For this transformation, we took z ¼ 4.52 nm (note that in the bulk, all positions are equivalent by definition). In addition, we calculated the free energy, DG z dp , for the same transformation at two different positions along the transport axis (z ¼ 1.02 nm near site Am1, and z ¼ À1.32 nm, near site Am5). This allowed us to shift the PMF profile for NH 3 with respect to that of NH 4 þ (which was previously shifted to have a value of zero in the bulk) before subtracting one curve from the other to obtain DG z dp (see Equation 1). Figure S3 shows two examples where the NH 3 PMF profile was shifted according to the free-energy calculations. Figure S3A shows the result acquired when statistics from only one channel were combined for the derivation of the NH 4 þ PMF profile, and Figure S3B shows the result for the NH 4 þ PMF from umbrella sampling the translocation across two different AmtB monomers (and combining the data from both monomers). The uncertainty in DG z dp at sites Am1 and Am5 determined by thermodynamic integration had an upper bound in uncertainty of 7 kJ/mol (see Table 1). At these positions, this translates to an upper bound of uncertainty in the resulting pKa of 1.2 units. In addition, we note that the NH 3 PMF profile could not be fit exactly to the points determined by thermodynamic integration (blue points with error bars in Figure S3). If, as an upper bound, we estimate the uncertainty to be 4 RT in the PMF profile, we obtain an additional uncertainty in the pK a profile of 1.7 units. Combining these estimates of uncertainty in the pK a calculation, we may give an upper bound in the determined pK a profile (of Figures 6B and S3C) of 2.9 units (or ;3 units, as mentioned in the text). The difference in the pK a obtained when utilizing 1-channel versus 2-channel sampling ( Figure S3) is obviously much greater than ;3 units in the hydrophobic lumen of the channel (;15-20 units around z ¼ 0.1 nm, see Figure S3C); however, in this region, the total pK a is about À45 units. Thus, regardless of the differences observed in this region, the apparent pK a of NH 4 þ is so low that it would be effectively impossible for Am to exist in its charged form. Despite any differences seen in the results from 1-channel and 2-channel sampling, we come to the same conclusion when determining the position where deprotonation and reprotonation of Am occur along its pathway through the channel (green vertical bars in Figures 6B and S3C).
All type-2 calculations were used to complete the occupancy state map in Figure 9A (free-energy values shown in blue). In these calculations, the NH 3 molecule (k ¼ 0) of interest at a particular site in the channel was transformed into a ''null'' molecule (k ¼ 1), where in the ''null'' state the NH 3 molecule had no interactions with any other molecule within the system and was harmonically restrained to its site. The value of h@H/ @ki P,T , at k ¼ 1, corresponding to the ''null'' state, was determined by linear extrapolation based on the tail-end values at k ¼ f0.85,0.90,0.95g. The free energy for each of these calculations is shown in Table 1. The free-energy values shown in the state map of Figure 9A were derived using the sequences of reactions shown in Figure S5B-S5E. Debye. The dipole moment of all system ions combined (i.e., the NH 4 Cl electrolyte) is seen to counteract the dipole of the AmtB trimer with an average value of hl z i ¼þ2600 6 1300 Debye. This value is precisely the same magnitude as the protein to within fluctuations due mostly to ion diffusion across the simulation cell boundary. (B) The system dipole moment. The value is zero to within fluctuations (hl z i ¼ 0 6 1500 Debye), due to the dielectric response of the electrolyte. Found at doi:10.1371/journal.pcbi.0030022.sg001 (2.2 MB EPS). Figure S2. Impact of AmtB versus POPC Membrane on System Ion Distributions Ion densities determined by calculating the average number of ions at a given position, z, that fall within 1.7 nm (in the xy-plane) of the center of mass of an AmtB monomer (black curves) and outside 3.7 nm (in the xy-plane) of the AmtB trimer center of mass (red curves). (A) Conditional NH 4 þ densities. The red curve shows normal behavior for NH 4 þ binding to the membrane surface. On the black curve, the region from z ; 1-2 nm, near the periplasmic vestibule, shows a concentration ;16 times larger than the region from z ; À1 to À2 nm, near the cytoplasmic vestibule. (B) Conditional Cl À densities. The red curve shows slightly less binding to the periplasmic membrane surface than the cytoplasmic surface. The black curve shows a concentration ;3 times smaller near the periplasmic vestibule than the cytoplasmic vestibule. (C) Illustration of the xy-planar regions for which the conditional density curves in (A) and (B) were calculated. The AmtB trimer is encircled by the gray cylinder with a radius of 3.7 nm whose center falls on the center of mass of the AmtB trimer. Three red cylinders of radius 1.7 nm encircle each AmtB monomer. Their centers coincide with the centers of mass of the three monomers. Found at doi:10.1371/journal.pcbi.0030022.sg002 (3.6 MB EPS). Figure S3. Analysis of the pK a Profile for NH 4 þ (A) PMF for NH 3 (red) and NH 4 þ (black) translocation as a function of the transport axis. The black curve is the result of umbrella sampling along z using a single monomer. The red curve was shifted with respect to the black curve based on thermodynamic integration calculations (datapoints shown in blue, see Table 1). A macroscopic view of the PMFs is shown in the inset. Note the large barrier for translocation (;150 RT) of NH 4 þ .

Supporting Information
(B) The same plot as shown in A, except the black curve is the result of sampling along z using two monomers (instead of only one). The major difference between the black curve here and that shown in A is the large barrier for NH 4 þ translocation (;100 RT, see inset). In either case the barrier represents an impassable region for NH 4 þ . Otherwise the curves are very similar. (C) The resulting apparent pK a for NH 4 þ as a function of the transport axis when using data from 1-channel statistics, as in (A) (dotted line), and 2-channel statistics, as in (B) (dotted line), to determine DG z dp in Equation 1. The vertical colored bars have the same meaning as in Figure 6. pK a ¼ 7 is marked with a horizontal dashed line. Note that despite the differences in the relative PMF shifts, the pK a profiles are nearly the same for all positive values, giving rise to precisely the same deprotonation region.  þ from our 25-ns production run, calculated as the average number of water molecules falling within the hydration shell boundary of NH 4 þ , at a given position. Hydration statistics were taken upon the condition that a NH 4 þ molecule was within 1.7 nm of the center of mass of a monomer backbone in the xy-plane (as was done in Figure S2 for partial-density calculations). Error bars were derived by dividing the trajectory into 500-ps blocks. Note that according to these data hydrated NH 4 þ does not pass the deprotonation regions at z ¼ 0.77 6 0.08 nm (periplasmic side) and z ¼À1.13 6 0.07 nm (cytoplasmic side) from the aqueous solution to enter the channel. The statistics for this curve were combined with those derived from umbrella sampling trajectories to produce the NH 4 þ hydration curve in Figure 6D.  The diagrammatic representation for a given occupied state is described in Figure 9. A filled circle within a site represents NH 3 occupancy, and an X within a site represents ''null'' occupancy (occupancy by a fictitious form of NH 3 restrained to the site, and having no interactions with the system). Free energies for each reaction (in kJ/mol) are written above the yield arrows, and are shown with their uncertainties from thermodynamic integration calculations in Table 1. The free energies resulting from the reactions (B-E) were used in the state map shown in Figure 9A (blue free-energy values). (A) The free energy to mutate NH 3 to a null molecule at site Am1. This value was used as the initial step for reactions (B-E), using the assumption that NH 3 at site Am1 is uncorrelated with the species in any luminal site. (B-E) Various reactions necessary for filling in the state map (blue free-energy values) of Figure 9A. Note that the free energy to move null NH 3 from one site to another is zero. All other free-energy values were derived from thermodynamic integration calculations (Table 1). Found at doi:10.1371/journal.pcbi.0030022.sg005 (672 KB EPS).