Modulating the Chemical Transport Properties of a Transmembrane Antiporter via Alternative Anion Flux.

ClC-ec1 is a prototype of the ClC antiporters, proteins that stoichiometrically exchange Cl- and H+ ions in opposite directions across a membrane. It has been shown that other polyatomic anions, such as NO3- and SCN-, can also be transported by ClC-ec1, but with partially or completely uncoupled proton flux. Herein, with the help of multiscale computer simulations in which the Grotthuss mechanism of proton transport (PT) is treated explicitly, we demonstrate how the chemical nature of these anions alters the coupling mechanism and qualitatively explain the shifts in the ion stoichiometry. Multidimensional free energy profiles for PT and the coupled changes in hydration are presented for NO3- and SCN-. The calculated proton conductances agree with experiment, showing reduced or abolished proton flux. Surprisingly, the proton affinity of the anion is less influential on the PT, while its size and interactions with the protein significantly alter hydration and shift its influence on PT from facilitating to inhibiting. We find that the hydration of the cavity below the anion is relatively fast, but connecting the water network past the steric hindrance of these polyatomic anions is quite slow. Hence, the most relevant coordinate to the PT free energy barrier is the water connectivity along the PT pathway, but importantly only in the presence of the excess proton, and this coordinate is significantly affected by the nature of the bound anion. This work again demonstrates how PT is intrinsically coupled with protein cavity hydration changes as well as influenced by the protein environment. It additionally suggests ways in which ion exchange can be modulated and exchange stoichiometries altered.

Introduction

The chloride channel (ClC) proteins constitute a family of transmembrane proteins in species ranging from prokaryotes to human beings that regulate a wide range of physiological processes.[1−5] Proteins within this family play a crucial role in the extreme acid resistance of bacteria (ClC-ec1)[5] and their genetic disruption or misregulation has been linked to multiple mammalian diseases, including myotonia congenita (ClC-1),[6,7] retinal degeneration (ClC-2),[4] Bartter’s syndrome type III (ClC-Kb),[2] Dent’s disease (ClC-5),[1] and osteopetrosis (ClC-7).[8,9] ClC proteins were initially assumed to be exclusively channels allowing the passive diffusion of Cl–, until ClC-ec1 was identified as an antiporter where Cl– and H+ ions are exchanged in opposite directions.[10] Importantly, this protein can also transport other anions, including NO3– and SCN–, with a similar conductance rate but altered proton coupling.[11−18] Because of their prevalence, transport diversity, and functional complexity, the ClC proteins have been aptly described as “a never ending source of surprises”.[19]

Among the members of ClC proteins, ClC-ec1 (Figure A) from Escherichia coli is the most structurally and mechanistically investigated.[10,20−24] It exchanges approximately two Cl– ions for one H+ actively and stoichiometrically under normal metabolic conditions.[10,24] As revealed by the wild-type (WT) crystal structure,[20,21] the cross-linked protein,[25] and the monomeric form,[26] ClC-ec1 is a homodimer in which two subunits function independently yet similarly. Though ClC-ec1 transports various anions, Miller and co-workers observed in electrophysiological experiments that the coupling of the polyatomic anions like NO3– and SCN– to H+ transport is weakened, or even extinguished, compared to Cl–/H+ coupling.[15] The stoichiometry of anion to proton flux is also shifted from (2.2 ± 0.1):1 for Cl–:H+ flow,[24] to 7–10:1 for NO3–:H+ flow, while proton transport (PT) is completely blocked by thiocyanate.[15] The goal of this work is to investigate these altered exchange processes in order to deepen our understanding of the ion exchange mechanism, the origin of chloride selectivity and the chemical facets that regulate exchange coupling and stoichiometry. The insights gained in this model system have broader implications for the driving forces of ion transport in other proteins, especially those that demonstrate partial or complete uncoupling effects, such as ClC-4 and ClC-5.[16,17,27]

Figure 1

Crystal structure of WT ClC-ec1 (PDB ID: 1OTS),[21] and two distinct states of Scen anion-binding for nitrate and thiocyanate. (A) The ClC-ec1 homodimer is shown with the key residues rendered and labeled in monomer A (blue). Approximate Cl– and H+ transport pathways are indicated by green and orange dashed arrows, respectively. Two bound chloride anions (green spheres) occupy the central (Scen) and internal (Sint) binding sites. The anion bound to Scen lies in the intersection of two ion pathways. (B) In the nitrate-bound protein, two states are characterized by the coordination of nitrate anion (van der Waals spheres) to S107 or Y445. (C) In the thiocyanate system, the state is differentiated by the orientation of the anion.

Unravelling the interaction between the transported ions in ClC-ec1 necessitates first understanding their binding sites and transport pathways. Bound to the crystal structure of ClC-ec1 are four Cl– anions, two in each monomer (see Figure A).[20] One of the two anions bound to each monomer is observed in the central binding site (Scen), located in the central cavity and surrounded by several conserved residues including S107, E148, I356, F357, A358, and Y445 (ClC-ec1 numbering). The Scen anion is stabilized by coordination to the backbone nitrogen atoms of I356 and F357, and to the oxygen atoms of S107/Y445 side chains. The other anion is bound to the internal site (Sint) with one side facing the internal bulk, and the other side coordinating the S107 and G108 backbone atoms. Crystal structures and simulations have suggested that the polyatomic anions occupy the same sites (Scen and Sint) during transport, but induce a significant shift in both local environment and proton transport functionality.[15,28] For the nitrate-bound ClC, this phenomenon is most likely attributable to the Scen site for which crystallographic and electrophysiological studies have demonstrated synergy with proton transport in ClC-ec1.[29] Similar synergistic interactions were also revealed in an eukaryotic Cl–/H+ exchanger.[30] Structural examination indicates that Scen is located at the intersection of the proton pathway and the anion pathway, additionally suggesting its potential effect on PT. This is further corroborated by the isothermal titration calorimetry (ITC) experiments that revealed the stability of Cl–, Br–, or NO3– bound to Scen in ClC-ec1.[31]

In contrast, the Sint-bound halide or nitrate observed in the internal part of the anion transport pathway has a low binding affinity.[31] This is not the case for SCN–, which is thought to have a higher affinity for the internal site based on undetectable electron density in Scen in a crystal structure of ClC-ec1 grown in SeCN– (a crystallographic analog of SCN–)[15] and ITC experiments.[31] We postulate that one of the reasons for uncoupling of H+ and SCN– transport could be the low affinity of SCN– to Scen, since the Scen-unbound ClC-ec1 was shown to be substantially less proton permeable in our previous publications.[32,33] However, that same work demonstrated how a Cl– anion at the Scen site has the largest impact on proton conductivity, and thus is central to Cl–/ H+ coupling. For these reasons, and because SCN– must pass through Scen with an unknown retention time, we focus herein on the impact of SCN– at Scen on PT. This also allows us to directly compare the effects of SCN– and NO3– to that of Cl– (all at Scen).

Though many biophysical studies of the anion transport mechanism have been performed,[29,31] the behavior of the transported proton has been more challenging to characterize experimentally. Site-directed mutagenesis experiments have revealed that two glutamate residues, E203 and E148, act as crucial proton loading sites for H+ transport.[34,35] However, between these two residues lies a ∼15 Å long hydrophobic cavity lacking in either ionizable residues or crystallographic water molecules, both of which are fundamental elements that enable PT. It therefore seems impossible that excessive protons can transport across such a dry area, when making judgements based solely on information from the crystal structure.

Despite the absence of a “proton carrier” in the crystal structure, the mystery in the hydrophobic cavity has been explained in a number of previous simulation studies in which the explicit PT is treated via a reactive molecular dynamics (MD) methodology that allows for Grotthuss proton shuttling between water molecules.[32,36−38] Moreover, nonreactive simulations have been used to predict that water molecules filled this hydrophobic region transiently depending on the protonation states of certain residues.[39,40] Multiscale reactive molecular dynamics (MS-RMD) simulations went further to demonstrate an important dynamical coupling of the PT with the increased hydration between E203 and E148.[32,37] This work also revealed that, in contrast to the expectation that PT from E203 to E148 should be rate-limiting and dependent on the presence of Cl– at Scen, the outward PT (E203 to E148) is actually relatively fast, even in the absence of Cl–cen.[32] Our subsequent work suggested that the role of Cl– bound to Scen during outward proton flux could instead be the facilitation of PT release from E148 to external bulk.[36] Our most recent work suggests that both of these PT steps may be rate-limiting and facilitated by Cl– occupancy in different contributing pathways to the total outward proton flux.[33] In contrast, the rate of proton uptake from external bulk to E148 was shown to be relatively fast at physiological pH.[33,36] This protonation does not require E148 to rotate up and is minimally affected by the presence or absence of an anion at Scen. The slowest rate for proton uptake to E148 was calculated to be (4 ± 1) × 102 ms–1 when Cl–cen is absent and E148 is rotated down,[36] which is more than 2 orders of magnitude faster than the Cl– turnover rate (2.3 ms–1) and PT rate (1.0 ms–1).[23,24] In the presence of Cl–cen, the proton uptake rate is slightly increased (less than 1 order of magnitude).[36] We further revealed that the inward proton flux is rate-limited by transport from E148 to E203, and dependent on the presence of Cl–cen.[33,36] Collectively, these previous studies explain the mechanism of ion exchange and the origin of the 2.2:1 ratio for Cl–:H+ flux. In the present paper, we consider the likely possibility that these rate-limiting steps are not transferable, and that PT within the central cavity may be significantly slowed by the presence of polyatomic anions at Scen. By slowing or blocking PT between E203 and E148, these anions could uncouple PT from anion transport since the requisite protonation of E148 for anion transport would still occur via protonation from the external side of the membrane. Therefore, we focus our attention in this work on PT between E203 and E148 as influenced by polyatomic anions.

Similar to the Cl–-bound system, the ease of hydrating the central hydrophobic region was investigated with standard MD simulations in the NO3–- and SCN–-bound systems; Jiang et al. observed shortened water networks within the central region of these two PT-suppressing/blocking systems.[28] However, analyzing the water network in the absence of an explicit hydrated excess proton (i.e., not including the delocalized protonic charge defect and Grotthuss proton shuttling) can present an incomplete, if not completely incorrect, picture of how hydration influences PT. The hydrated excess proton has been clearly shown to significantly alter water structure and dynamics,[32,41,42] including within ClC-ec1.[32] Thus, the causality between the change in water structure and the decrease in PT rates warrants further exploration, by explicitly including the dynamic coupled behavior of proton transport and water hydration in these polyatomic anion-bound systems.

The investigation of the explicit proton transport as affected by polyatomic anions is challenging in three distinct ways: (1) as noted earlier the explicit excess proton moves along the water chain through consecutive O–H bond breaking and reforming steps (Grotthuss mechanism), which are not handled correctly in standard nonreactive MD simulations; (2) as was proposed by Jiang et al.,[28] these polyatomic anions in the central cavity, unlike monatomic Cl–, show two distinct binding modes (Figure B,C), which must be treated carefully since the environment within the cavity will clearly influence the PT rate;[32] and (3) the polyatomic anions might even interact with excess protons directly in the protein environment during PT between E203 and E148.

In this work, we demonstrate how ion exchange is modulated by the chemical nature of the anion and the protein environment. We were able to overcome the challenges above by using MS-RMD simulations and electronic structure calculations to investigate the influence of polyatomic anions on PT. Quantitative free energy calculations of explicit PT and the associated changes in hydration reveal a novel three-phase mechanistic coupling between hydration and PT in the central cavity. As opposed to simply the frequency of water filling the protein cavity, we find that the hydrogen-bond connectivity in narrow regions of the cavity is the critical feature that influences PT and ion coupling. Thus, PT through the central cavity can be either facilitated (Cl–) or inhibited (polyatomic anions) depending on the anion’s size, proton affinity, protein interactions, and influence on hydration. Surprisingly, we find the proton affinity of the anion (i.e., the pKa of its conjugate acid) has little influence, while its size, stability and influence on hydration are the dominant forces of change. These findings help to explain the coupling mechanism in ClC-ec1, and more broadly suggest ways in which ion exchange can be modified.

Results and discussion

Anion Basicity and Binding Modes

We first determined whether or not the polyatomic anions would interact directly with the proton (i.e., be protonated) as it is transported through the central cavity. Extensive hybrid quantum mechanics/molecular mechanics (QM/MM) MD simulations were performed and ruled out such possibility, which is described in the Supporting Information (SI). Since we observed that the anions stayed stable in Scen for ∼0.5 μs of MD simulations, and that they bound to the protein in two different states (Figure B,C), we next wanted to assess whether or not the two anion binding modes needed to be explicitly sampled in our PT analysis. Thus, we analyzed the relative stability and rate of transitioning between the two modes with MS-RMD umbrella sampling simulations including the explicit H+ on E203. The potentials of mean force (PMFs, i.e., free energy profiles) for the motion of anions (Figure ) confirmed two distinct binding modes for both anions and revealed a larger transition barrier for thiocyanate. On the basis of Figure B,C, we defined state I as the left minima and state II as the right. In the nitrate case, the anion either coordinated with S107/Y445 (state I) or not (state II). For the SCN–-bound protein, either the S-side (state I) or the N-side (state II) of the anion interacted with G149/F357, while the other side coordinated the lower part of the cavity (S107/Y445). Contrary to previous results reported,[28] we found that state II for both nitrate and thiocyanate was thermodynamically favorable. This difference is likely due to the different force fields and methods employed. Jiang et al.[28] used a nonreactive TIP3P water model with no enhanced sampling, whereas we performed reactive MD with an excess proton on E203 and the SPC/Fw water model[43] as well as umbrella sampling (see Methods). Given the sensitivity of water dynamics to the underlying model, the water environment around the anion, and hence the relative enthalpic and entropic stability, could also be affected by the water force fields employed. Observing that the nitrate should experience less steric hindrance after the S107 side chain rotates down from state I to state II, we suggest that the unexpected relative stability of nitrate state II is due to entropic freedom in a larger free space, compensating for the energetic loss of coordination with the S107 and Y445 side chains. The free energy barrier and rate constants (Table S2) extracted from the PMFs indicate that the rotation of nitrate is relatively fast (picosecond time scale) and able to be sampled within the production runs (nanosecond time scale), whereas the flipping of the bulkier thiocyanate anion is slower (1.1 ns) and will not be sampled sufficiently. The fast state transition of NO3– is due to its relatively small size as compared to the protein cavity, and the relative weakness of both steric and coordination effects. In comparison, the slow state transition of SCN– is due to its linear and more bulky structure. The surrounding residues squeeze SCN–, disfavoring intermediate configurations and slowing the flipping motion of SCN–. On the basis of these results, we decided to treat the two thiocyanate binding modes separately in the PT analysis, whereas the faster transitions for nitrate were deemed sufficient to be captured in PT simulations without distinguishing or explicitly sampling the two binding modes.

Figure 2

PMFs with error bars for the two binding modes of nitrate (A) and thiocyanate (B), within the central cavity of ClC-ec1. “TS” is short for transition state. Screenshots of the representative configurations of the metastable states, i.e., local minima, are depicted. The collective variable (CV) for the PMF is defined as (A) the smallest distance between the nitrate oxygens and the oxygen on the side chain of S107 residue; (B) dS,ref – dN,ref, the difference in the two distances from “ref” to either “S” or “N”, in which the subscripts “S” and “N” stand for the sulfur atom and the nitrogen atom of the thiocyanate anion in the PT pathway, respectively, and the “ref” represents the midpoint of the α-carbons of the residues R147 and A358.

Proton Transport Simulations

Our initial studies described above determined the need for three simulations of PT: two for SCN– and one for nitrate. Preliminary 1D sampling of the PT across the central cavity showed a significant drop of the hydration level at the state of protonated glutamic acid, implying the coupling between PT and the water environment. Hence, 2D umbrella sampling simulations were performed to construct PT PMFs. The PMFs and calculated rate constants are summarized in Figure and Table , respectively. The increased free energy barrier and calculated rate constants agree with the experimental results that PT is partly uncoupled to the nitrate transport and completed blocked by the thiocyanate.[15] Combining with the experimental measurement where the ratio of anion transport rates of Cl–:SCN–:NO3– is approximately 3:3:1[15] and the Cl– turnover rate of 2.3 ms–1,[23,24] the shift in PT rate is qualitatively consistent with the stoichiometry change from 2.2:1 for Cl–:H+ flow,[24] to 7–10:1 for NO3–:H+, and to almost infinity for SCN–:H+.[15] The consistency justifies our computational results and allows for additional analysis based on these 2D PMFs.

Figure 3

2D PT PMFs for nitrate (A) and thiocyanate in binding mode I (B) and II (C) tracking the improved definition for the ratio-based and unitless CEC position (horizontal axis) as defined in eq and water density (vertical axis) defined in the SI. Irrelevant high-energy areas (red) are not sampled to reduce the cost of computation. The black lines trace the minimum free energy paths for PT from E203 to E148. (D, E, and F) The extracted 1D PMFs, along the minimum free energy path, from the corresponding 2D PMFs above. The error of the PMFs is ∼0.2 kcal/mol for the nitrate case and ∼1 kcal/mol for the thiocyanate case, estimated from splitting the trajectories and block-averaging. Key points are labeled in each PMF denoting: (A) protonated E203; (B) the proton in a hydrated channel just above E203; (C) the transition state; (D) the proton passing through a connected water network to E148; and (E) protonated E148 and dehydration of the channel. Note the much larger energy scale of the vertical axis in panels D, E, and F versus that of the same lettered panels of Figure S1 for the WT protein.

Table 1

Summary of PT PMFs with Different Anionsa

System	ΔF‡ (kcal/mol)	k (s–1)	kref (s–1)
Cl–	7.4 ± 0.3	(7 ± 3) × 107	4 × 108
NO3–	14.6 ± 0.2	14 ± 5	∼7 × 101
SCN–, I	21 ± 1	(2 ± 3) × 10–4	<1
SCN–, II	28 ± 1	(1 ± 2) × 10–7	<1

ΔF‡ is the energy barrier for the proton transport between E203 and E148, from the higher-energy state (E203 for Cl–/ NO3– and E148 for SCN–) to the lower-energy state (E148 for Cl–/ NO3– and E203 for SCN–), and k stands for the rate constant calculated from transition state theory. The experimentally measured anion rate constant combined with the anion/proton stoichiometry is used to calculate the reference kref for NO3–, and kref for SCN– is based on the detection limit of ion transport.[15,23,24] Since PT in the central region is not rate-limiting in ClC-ec1, our previous result for PT from E203 to E148 is listed as the reference kref for Cl–.[32]

To understand the PT mechanism, 2D PMFs as well as 1D traces along the minimum free energy paths are inspected. According to the 2D PMFs, the PT process can be categorized into three phases: (1) hydration (path A to B in PMFs), (2) proton transport (path B to D in PMFs), and (3) dehydration (path D to E in PMFs). In the hydration phase, the hydrophobic cavity is solvated by 3–5 water molecules, as shown by the increase in water density collective variable (CV) along the initial part of the minimum free energy paths. For both nitrate and thiocyanate, the initial deprotonation of E203 is facile once the hydration level reaches 3–3.5. However, for nitrate the cavity can easily fill to ∼5 waters before or during E203 deprotonation. In contrast, it is harder to increase the solvation to ∼5 or more water molecules in the thiocyanate system, requiring more free energy and coupling to PT (meaning the excess proton transitioning to the water above E203 facilitates more water filling the cavity). In the proton transport phase, the excess proton progresses toward E148. This is an uphill process in free energy, peaking as the proton passes the anion in each case. For thiocyanate, it is coupled to increasing hydration. Finally, E148 is protonated and the cavity returns to the low hydration state due to the absence of an excess proton. These results are consistent with our preliminary 1D sampling.

Based on several prior computational results, it was expected that the hydration phase should contribute to the proton transport free energy barrier most significantly,[28] and that the proton transport should couple with the hydration.[32,41] Surprisingly, these two CVs do not couple clearly in terms of the minimum free energy path, especially for the nitrate system (Figure A), and the hydration phase does not induce much of a barrier (less than 5 kcal/mol in the nitrate case). This is attributable to the fact that the water density CV ξ2 introduced in our previous paper[41] simply reflects the water count of a specified region, but not necessarily the water wire connectivity with the excess proton in it. In the hydration phase, water molecules cluster below the anion, which is relatively easy due to the proximal ionizable residue (E203). Importantly, cavity hydration is not necessarily always equal to water wire connectivity.

In order to better understand the origin of the dominant contribution to the PT barrier, we analyzed the correlation between free energy and water connectivity, by plotting the extracted 1D PMF along with either the probability of forming a continuous water-wire (Figure A,B,C), or the ratio of the length of the disconnected area to the distance between two glutamates, E203 and E148 (Figure D,E,F). These plots indicate a strong correlation between the free energy barrier and the formation of a continuous water-wire, as is highlighted by two facts: (1) the water network is only connected around the dividing surface where the derivative of the free energy over the path is near zero; (2) the free energy and the ratio of disconnected distance to total distance is inversely related. This explains the unexpected free energy barrier in the PT phase and recaptures the correlation between PT and water dynamics. Thus, it is the water connectivity along the PT pathway in the presence of the explicit excess proton that is most relevant to the PT barrier. Interestingly, the probability of forming a continuous water-wire in the state II of thiocyanate system (Figure C) is extremely low, indicating that the water-wire actually breaks on one side when it is forming on the other side. According to the correlation between water-wire formation and PT barrier, this result is consistent with the highest PT free energy barrier in the system.

Figure 4

Overlap of extracted 1D PMF (black) and either the probability of forming a continuous water-wire (A, B and C) or the ratio of the length of the disconnected area to the distance between two glutamates, E203 and E148 (D, E and F) (blue), plotted against the minimum free energy pathway. They are plotted for the nitrate-bound protein (A, D), or the thiocyanate-bound protein in state I (B, E) or state II (C, F).

Water-Network Analysis

For each of the polyatomic anion PMFs, the free energy for PT peaks (at the transition state) as the proton passes the bulky anion. This was not the case for Cl– (see Figure S1). At the initial stage of the PT phase where free energy escalates, the water network is incomplete (Figure A,D). Connecting the water network to form a so-called “water-wire” past the bulky anions is strongly coupled to, and enabling of, PT past the anions to E148. Hence, the free energy ceases to increase when the water network around the bound anion is connected (Figure B,E). The position of center of excess charge (CEC, defined in eq ) is exactly next to the anion at this stage. Finally, the system is stabilized after the proton passes the anion, binds to E148, and the water-wire is again disconnected (Figure C,F). For the entire PT process, only when the excess proton is near the anion does a continuous water-network (i.e., continuous near the central anion) become significantly observable, which is near the dividing surface. Similar proton-induced wetting has been observed in other biological systems and materials.[41,42,44] The distinct state change of the water-wire near the dividing surface confirms that the water molecules close to the polyatomic anion are energetically unfavorable. Thus, the coupled phenomena of PT and hydration, as well as the challenge of moving a proton through a confined space with limited delocalization, significantly increases the free energy barriers for PT between E203 and E148. Based on this analysis, the dynamic process of the PT phase can be more accurately summarized as (1) proton transport from E203 to below the anion along with the elongation of the water-network but no connectivity past the anion (state B); (2) the water network connecting past the polyatomic anion, in the presence of a proximal excess proton, enabling PT past the anion (state C); and (3) breaking of the water network with the loss of the excess proton as it binds E148 (state D). The dividing surface lies in the second stage.

Figure 5

Representative configurations of the nitrate-bound (A, B and C) and thiocyanate-bound (D, E and F) systems in the initial stage of the PT phase (A, D), connected water network (B, E), and final stage of the PT phase (C, F). The Scen anion is shown in VDW representation, and the hydronium/protonated glutamic acid is rendered orange for clarity. Key residues (labeled in panel A) and water molecules are colored according to element.

To further elucidate the molecular basis of the PT/water-induced barrier, we analyzed the anion structure close to the free energy barrier and found significant steric hindrance by the polyatomic anions. Unlike monatomic Cl–, polyatomic anions like NO3– and SCN– occupy more space within the central region. NO3– has a trigonal planar geometry, with the negative charge distributed over three O atoms. The effective radius of NO3– is ∼1.97 Å, which is larger than that of Cl– (1.64–1.82 Å).[45,46] A slightly larger spatial occupation of nitrate anion is enough to slow PT across the central region, but insufficient to completely block the H+ transport. In contrast, SCN– has a linear structure that spans 4.77 Å in length,[47] making it much bulkier than either Cl– or NO3–. The increased size of SCN– increases the free energy necessary to connect the water network and transport the excess proton through the transient water chain. The influence of steric hindrance is consistent with our finding that the spatial effect slightly outcompetes the coordination effect in the thiocyanate motion PMF (Figure ).

Conclusions

We have performed extensive multiscale reactive MD simulations combined with free energy sampling to investigate the influence of polyatomic anions on PT and the related mechanisms of ion coupling and stoichiometric ion exchange in ClC-ec1. In our simulations, PT was treated explicitly, including Grotthuss shuttling and proton charge defect delocalization. We first quantified the relative stability and transition rates between two binding modes of the nitrate and thiocyanate anions in the central binding site of ClC-ec1, Scen. We found that the anions are stable for at least the submicrosecond time scale in the central site. Our PMF calculations reveal that separate treatment of the two thiocyanate binding modes was necessary due to the slow interconversion. In contrast, nitrate rotates relatively quickly, making it possible to sufficiently sample this transition during a single PT PMF calculation. This difference between nitrate and thiocyanate can be explained by the increased steric restriction of the bulky thiocyanate. We then defined and justified a new reaction coordinate to characterize PT between E203 and E148, which incorporates the excess proton CEC and the distance between the two glutamate residues simultaneously (see SI). This CV significantly reduces the computational demands. It not only captures two crucial factors of PT in a single variable, but also introduces no new spurious effects on the PMF based on comparisons to previously published PMFs[32] for PT for Cl–-bound ClC-ec1. Next, we used QM/MM simulations to verify that the polyatomic anions do not react with proton directly during the PT process in the protein, and hence do not require explicit MS-RMD treatment (SI).

After addressing the challenges of studying the NO3–/SCN–-bound ClC-ec1, we constructed three 2D PMFs for the PT across the central region, one for NO3– and two for SCN–. The rate constants calculated from the PMFs utilizing transition state theory are consistent with the experimental measurements, in which PT is partly uncoupled to the nitrate transport and completely uncoupled to thiocyanate transport. A previously published hypothesis about the water structure being influenced by the polyatomic anions was confirmed but significantly expanded to a three-phase mechanism of proton shuttling, on the basis of quantitative PT free energy calculations. Interestingly, these PMFs revealed that simply increasing the number of water molecules in the region of PT is not sufficient to enable PT. The critical hydration CV seems to be water network connectivity, which we demonstrate by analyzing the probability of water connectivity and the ratio of the length of the disconnected area to the distance between two glutamates. Both of these measures are highly coupled with PT through the transition state, as indicated by the correlation with the 1D PMFs for all three systems (one for NO3– and two for SCN–). A continuous water-wire with an explicit excess proton in it is in fact able to exist in the NO3–/SCN–-bound protein, albeit transiently and conditionally, contrary to the hypothesis arising from a prior standard empirical MD simulation in which the explicit PT was not treated.[28]

In-depth analysis on the PT phase reveals three stages for the process. First, there is an increase in hydration weakly coupled to PT to just below the anion binding site. Second, the transition state is defined by strongly coupled connection of the water network and PT past the anion. The cause of the free energy barrier, and thus decreased proton coupling, is the steric hindrance and unfavorable hydrogen bond connectivity induced by NO3– and SCN–. For the smaller polyatomic anion NO3–, the PT free energy barrier is raised enough to make PT through central cavity clearly rate-limiting and slower than it is with Cl–, whereas the bulkier SCN– completely inhibits PT. Third, the water network breaks as E148 is protonated. This mechanism demonstrates that it is not primarily the hydration of the cavity that is rate limiting, but rather forming a connected hydrogen bonded pathway past the steric hindrance of the polyatomic anions. This provides insightful contrast to the common assumption that blocking water influx pathways is an effective way to inhibit PT (or an explanation of inhibited PT). For this system, and we anticipate for many others, water influx is a relatively fast and easy process compared to water connectivity through narrow regions. The latter was in fact recently found to be the cause of the rate limiting step for PT in cytochrome c oxidase.[42,48]

Collectively, this work shows how seemingly subtle changes in the chemical nature of the anion can flip it from facilitating coupled proton exchange, to inhibiting or even blocking it in both directions. The shift in coupling changes the associated ion exchange stoichiometries, instead of completely blocking anion transport, because the rate and ease (i.e., pKa) of protonating E148, which is essential for anion flux, from the external bulk is fast and minimally influenced by the presence of the anions.[33,36] Rather, PT is hindered/blocked in the central cavity between E148 and E203. The demonstrated entanglement between PT, water dynamics, and the local environment is expected to play a significant role in a number of other proton-conducting proteins. Our insights also potentially suggest guidelines for one way of controlling the function of the ClC family of proteins by influencing hydration, either through tuning the spatial constriction of the central cavity via binding, or by adjusting the steric hindrance caused by the channel-forming residues. However, our conclusions for these permeant anions are not necessarily transferrable to the uncoupling of PT to fluoride-bound ClC, due to the impermeant nature[3,49] and the strong electronegativity of F–. Therefore, additional research on PT in fluoride-bound ClC proteins, including a ClCF exporter,[50] is needed and will be the focus of future efforts. Additionally, we hope to identify new hydration CVs that directly track water connectivity and to investigate the Sint-bound SCN– system, as well as the complete anion/proton transport cycle employing multiscale kinetic modeling (MKM).[33]

Methods

The employed method for simulating explicit proton transport, MS-RMD, performs an on-the-fly diagonalization of a Hamiltonian matrix to identify a linear combination of “states”, each with a different bonding topology, that delocalizes the excess proton’s positive charge and enables Grotthuss shuttling without running much more expensive QM/MM MD simulations. The procedures for setting up the system, a more detailed introduction to MS-RMD, and the parametrization procedure of MS-RMD models[51−54] for residues E148 and E203 are explained in detail in the SI. Initial configurations for the simulations were mutated from a structure obtained from a previous study of this system,[32] which was equilibrated from the ClC-ec1 crystal structure (PDB ID: 1OTS).[21] After equilibration for ∼0.5 μs, the production MS-RMD umbrella sampling runs were performed with the RAPTOR software[51] embedded in the LAMMPS MD package (http://lammps.sandia.gov),[55] with umbrella potentials implemented in the PLUMED package.[56] Details of generating PMFs are described in SI for the sake of simplicity.

Simulations of the Anion Motions

The relative stability of the central anions in each binding mode was investigated by calculating the PMF for the rotation or flipping of the anion as described by a single continuous CV. In the nitrate-bound system, the CV definition was inspired by the hydrogen-bond interaction between the residues and the anion.[28] Among the candidates, the selected CV, smallest distance between the nitrate oxygens and the oxygen on the side chain of S107 residue, shows distinct bimodal distribution within the classical simulations. For the thiocyanate-bound system, we first defined a reference point located slightly above the anion. The midpoint of the α-carbons of residues R147 and A358 was chosen for stability. Then the anion’s flipping motion was characterized by the difference in two distances, from the reference point to either sulfur atom or the nitrogen atom of the thiocyanate at Scen. With the CVs defined above, umbrella sampling for ∼10 ns was performed to construct the PMFs.

Simulations of the Proton Transport

Four 2D PMFs (two for the SCN– cases, one for NO3–, and one for WT) were calculated from umbrella sampling simulations. Each of them is described with two CVs, one (ξ1) representing the extent to which excess proton CEC transports, and the other one (ξ2) indicating the presence of water molecules surrounding the transport channel.[32]

The first CV is the ratio of the distance between the CEC and the starting point to the distance between the end points. The advantage of this CV is that it quantitatively captures the movement of the CEC and the relative distances between donor and acceptor simultaneously. Mathematically, the first CV is described aswhere rCEC is a prescribed unit vector pointing from E203 to E148, rE203, and rE148, are the coordinates of the CEC, i-th oxygen of E203 side chain, and j-th oxygen of E148 side chain, respectively, with i and j enumerating from 1 to 2. The min function used in ξ1 is a smoothed approximation (defined in SI) of the mathematical minimum function. This CV, different from the curvilinear path CV used in our previous work,[32] describes the system in a more systematic and complete manner. See SI for more discussions.

The second CV (ξ2) is the effective number of water molecules in a box, with the center of the box and the principal axis defined through our published procedure,[32] and further described in the SI. The boundary of the box was determined so that waters outside the central region are not counted.

Rate Constant Calculations

As is derived from the transition state theory,[57,58] the rate constant k was calculated fromin which T represents the system temperature, and the effective mass meff is estimated from meffv2 = kBT based on the equipartition theorem where the instantaneous velocity v of the excess proton CEC motion is calculated from the finite-difference method. In the equation above, V(x) stands for system potential energy, and x is the full set of system variables. The probability (Boltzmann factor) is integrated over the dividing surface S in the numerator and over the reactive basin A in the denominator. The characteristic time was calculated as the inverse of rate constant.

Water-Wire Analysis

A continuous water-wire is defined as a water network starting from E203 to E148, connected through hydrogen bonds, determined by two criteria: (1) the distance between the donor and acceptor heavy atoms is shorter than 3.0 Å; and (2) the angle formed by the donor atom, the central hydrogen, and the acceptor atom is larger than 150°. The breadth-first algorithm was used to search for any possible continuous water network. In the case of nonconnecting water-wire, the searching algorithm was performed twice, one from E203 to E148 and the other in the opposite direction, with maximum depth set as 20. The smallest distance between two groups of heavy atoms in the water networks starting from either E203 or E148 is divided by the smallest distance between two groups of side-chain oxygen atoms of either E203 or E148, and plotted against the minimum free energy path.