Molecular Mechanism of ATP Hydrolysis in an ABC Transporter.

Hydrolysis of nucleoside triphosphate (NTP) plays a key role for the function of many biomolecular systems. However, the chemistry of the catalytic reaction in terms of an atomic-level understanding of the structural, dynamic, and free energy changes associated with it often remains unknown. Here, we report the molecular mechanism of adenosine triphosphate (ATP) hydrolysis in the ATP-binding cassette (ABC) transporter BtuCD-F. Free energy profiles obtained from hybrid quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations show that the hydrolysis reaction proceeds in a stepwise manner. First, nucleophilic attack of an activated lytic water molecule at the ATP γ-phosphate yields ADP + HPO4 2- as intermediate product. A conserved glutamate that is located very close to the γ-phosphate transiently accepts a proton and thus acts as catalytic base. In the second step, the proton is transferred back from the catalytic base to the γ-phosphate, yielding ADP + H2PO4 -. These two chemical reaction steps are followed by rearrangements of the hydrogen bond network and the coordination of the Mg2+ ion. The rate constant estimated from the computed free energy barriers is in very good agreement with experiments. The overall free energy change of the reaction is close to zero, suggesting that phosphate bond cleavage itself does not provide a power stroke for conformational changes. Instead, ATP binding is essential for tight dimerization of the nucleotide-binding domains and the transition of the transmembrane domains from inward- to outward-facing, whereas ATP hydrolysis resets the conformational cycle. The mechanism is likely relevant for all ABC transporters and might have implications also for other NTPases, as many residues involved in nucleotide binding and hydrolysis are strictly conserved.

Introduction

Nucleoside triphosphate (NTP) hydrolysis is one of the most fundamental chemical reactions in biology.[1,2] The free energy release upon cleavage of NTP into nucleoside diphosphate (NDP) and phosphate (P) is used by various proteins such as motors,[3,4] signaling proteins,[5] and transporters.[6] The basis for the function of these proteins is their high catalytic activity: some NTPase enzymes accelerate hydrolysis by up to 9 orders of magnitude over the reaction in aqueous solution. However, despite its importance, a detailed atomic-level understanding of the molecular mechanisms that underlie phosphoryl transfer reactions is still partly lacking.[1,2] This discrepancy can at least partly be attributed to experimental challenges linked to achieving the high spatial and temporal resolution required for studying this intricate chemical reaction in biomolecular systems.

The hydrolysis reaction involves several key steps, including nucleophilic attack of a lytic water molecule at the γ-phosphate, cleavage of the scissile phosphate bond, proton transfer (PT) events, and hydrogen bond rearrangements. Theoretical methods such as hybrid quantum mechanical/molecular mechanical (QM/MM)[7] and ab initio molecular dynamics[8] approaches can provide the lacking atomic-level insights into the structural and energetic details of chemical reactions in complex environments. Hence, such methods have been applied to study phosphate hydrolysis in aqueous solution,[9−15] guanosine triphosphate (GTP) hydrolysis in signal transduction proteins and translational GTPases (recently reviewed by Carvalho et al.[16]), and adenosine triphosphate (ATP) hydrolysis in molecular motors such as kinesin,[17] myosin (recently reviewed by Kiani and Fischer[18]), actin,[19,20] and F1-ATPase.[21−24] These studies focused on the following questions. First, does hydrolysis proceed via an associative, dissociative, or concerted mechanism? These pathways are (i) idealized, and (ii) energetically very close, which has led to controversial discussions in the literature.[2,15,16,25,26] Second, what is the nature of the proton-accepting group, and along which pathway(s) does proton transfer occur? A third question concerns the height of the activation energy barrier, and whether the hydrolysis reaction inside the protein is ex- or endergonic. The free energy of ATP hydrolysis is ΔG = −7.3 kcal/mol in aqueous solution under standard conditions, and even more exergonic under physiological conditions.[27] However, ΔG of the hydrolysis reaction itself might be different in the protein environment. For example, ΔG is close to zero in F1-ATPase.[28]

None of the above questions have been comprehensively answered for ATP-binding cassette (ABC) transporters, which couple the free energy gain of ATP binding and hydrolysis to the transport of substrate molecules across biological membranes.[6,29−32] All ABC transporters share the same architecture, a dimer of two nucleotide-binding domains (NBDs) or ATP-binding cassettes (ABCs) that bind and hydrolyze ATP, and two transmembrane domains (TMDs) that provide the pathway for substrate translocation across the membrane. The NBDs are usually highly conserved among different ABC transporters, suggesting a common ATP hydrolysis mechanism. ABC transporters can be understood as ATP-powered molecular machines. Their working cycle is based on the conversion of the chemical energy stored in ATP into mechanical work, which needs to be transmitted from the NBDs to the TMDs to induce the conformational changes necessary to translocate transport substrates. Whether the actual “power stroke”, in terms of the release of free energy that is necessary to drive the conformational changes, is provided by hydrolysis or binding of ATP remains an open question, the answer to which might depend on the type of ABC transporter.[33−35] X-ray crystal structures, including some with transition state analogs, yielded valuable insights into possible mechanistic scenarios of ATP hydrolysis in ABC transporters.[36] However, X-ray crystallography can only provide static snapshots and cannot resolve the positions of the protons. In addition, the energetics of the reaction, in terms of free energy profiles of the individual reaction steps, remains unknown. The link between several human diseases and the malfunction of ABC transporters[6] further underlines the significance of understanding their mechanisms of action.

To our knowledge, ATP hydrolysis in ABC transporters has been addressed in only three theoretical studies published so far. Figure S1 shows a schematic overview of the different mechanistic scenarios proposed. In the first study, Zhou and co-workers[37] used the AM1 semiempirical quantum chemical method within a QM/MM framework to investigate ATP hydrolysis in the Escherichia coli hemolysin toxin transporter HlyB. On the basis of potential energy surface calculations, they concluded that the conserved H-loop histidine residue (H191 in Figure , the corresponding residue in HlyB is H622) is directly involved in the hydrolysis. According to Zhou and co-workers, the H-loop histidine first acts as a general acid by donating its NϵH proton to the γ-phosphate, and subsequently accepts a proton from the lytic water molecule.[37] This proposed mechanism is surprising, as it proceeds via a negatively charged histidine imidazole ring, which at pH 7 is unfavorable in free energy. Indeed, the activation energy barrier of 22 kcal/mol obtained from the calculations[37] is rather high. In addition, the highly conserved Walker B glutamate (E159 in Figure ), mutations of which impair ATPase activity and which is thought to be the catalytic base,[36,38−40] is not directly involved as a general base in this proposed mechanism. The other two studies,[41,42] which investigated the maltose transporter MalK, also suggest a direct involvement of the H-loop histidine (H192 in MalK). Huang and Liao employed a QM/MM nudged elastic band method to calculate energy profiles, using the B3LYP density functional for the QM region.[41] Assuming a positively charged (doubly protonated) H-loop histidine, they propose that two protons are transferred upon nucleophilic water attack at the γ-phosphate, one from H192 to the γ-phosphate and the second from the lytic water to the conserved glutamate. These two proton transfers occur at the same point along the chosen reaction coordinate, and the reaction proceeds via a single barrier of 19 kcal/mol. These findings were confirmed by Hsu and co-workers,[42] who carried out QM/MM metadynamics simulations using the BLYP density functional; however, the reported reaction barrier of 10.5 kcal/mol is rather low. Furthermore, in this mechanism the protonation states of the H-loop histidine and the conserved glutamate at the end of the reaction differ from the beginning, and resetting the catalytic cycle would thus require additional protonation/deprotonation events. Another peculiarity of this second proposed mechanism is that a doubly protonated H-loop histidine might only be populated at pH < 7. However, at such low pH, ATP could be protonated as well,[16] which would have implications for the mechanism. In summary, the mechanism of ATP hydrolysis in ABC transporters remains elusive.

Figure 1

ATP hydrolysis mechanism in the ABC transporter BtuCD-F. (A–C) Snapshots from QM/MM MD simulations of the reactant state (ES), transition state (TS1), and product state (IS1) of the first reaction step. (D) Potential of mean force (PMF) along the Pγ–OW distance from the simulations with the standard and enlarged QM subsystems (black and gray lines, respectively); the larger QM subsystem additionally includes the side-chains of H191 and K39. (E) Distance between the transferred proton of the attacking water molecule, HW, and the Oϵ atom of E159. (F) Change of the breaking Pγ–Oβ bond distance (upper panel) and distance between the proton-accepting Oϵ atom of E159 and the transferred water proton (lower panel) along the reaction coordinate. The black and gray lines were obtained with the standard and enlarged QM subsystem, respectively. (G) Projection of the ES (black), TS1 (magenta), and IS1 (cyan) trajectories onto the Pγ–Oβ and Oϵ–OW distances.

Here, we use QM/MM molecular dynamics (MD) simulations to investigate the structural and energetic details of the ATP hydrolysis reaction in the vitamin B12 importer BtuCD-F from E. coli, one of the structurally most thoroughly characterized ABC transporters.[43−47] Unlike most previous QM/MM studies, the ATP hydrolysis mechanism was not investigated with minimum energy path calculations, but on the basis of potentials of mean force (PMFs) at room temperature. The MD sampling necessary for obtaining statistically meaningful PMFs makes this approach computationally much more expensive (the accumulated sampling time of our free energy simulations is more than 1 ns, which with the 1 fs time step used to integrate the equations of motion in the MD corresponds to more than 1 million electronic structure calculations), but it takes entropic contributions explicitly into account. Previous studies have shown that entropy can play a significant role for the free energy profile of NTP hydrolysis reactions.[18,48] In the QM/MM MD simulations, the MPW1B95 density functional[49] with the 6-31+G(d,p) basis set was used for the QM part, which consists of the triphosphate chain, Mg2+ ion, and several protein side-chains and water molecules (see the Methods section). MPW1B95 was chosen, because it can very accurately describe the energetics of the hydrolysis reaction, as shown in a benchmark study of phosphodiester hydrolysis that compared 52 density functionals.[50] When used together with a basis set that includes diffuse functions on heavy atoms, MPW1B95 energies deviate by less than 2 kcal/mol from approximated CCSD(T)/CBS values.[50]

Our results show that ATP hydrolysis in the ABC transporter follows a three-step mechanism. In the first step, nucleophilic attack of the lytic water molecule at the γ-phosphate cleaves the phosphate bond and yields ADP + HPO42– as intermediate product. Upon formation of the new Pγ–O bond, the attacking water molecule transfers a proton to the highly conserved glutamate E159, which acts as a catalytic base. In the second reaction step, the proton is transferred back from E159 to HPO42–, leading to ADP + H2PO4–. Finally, in the third step, the hydrogen bond network rearranges such that H2PO4– is stabilized by hydrogen bonds to ADP and E159; this rearrangement also involves changes of the Mg2+ ion coordination. From our computed free energy profiles, we estimate a reaction rate constant of ca. 7 s–1, in good agreement with the experimental value of ca. 0.5 s–1. The Δ of the overall reaction is close to zero, ΔG = +1.8 kcal/mol. Hence, cleavage of the phosphate bond itself cannot provide the power stroke for conformational changes of the transporter that are linked to substrate transport. Rather, the free energy of ATP binding is essential for inducing the tight dimerization of the NBDs and the associated transition of the TMDs from an inward- to an outward-facing conformation, whereas ATP hydrolysis is required to reset the conformational cycle.[51−53] Furthermore, although not explicitly investigated in this work, we speculate that unbinding of the hydrolysis products P and ADP likely plays an important role as well.

Results and Discussion

Attack of the Lytic Water Molecule

The first step of the hydrolysis reaction is the nucleophilic attack of the lytic water molecule at the γ-phosphate of ATP and breakage of the scissile Pγ–Oβ bond. As a prerequisite, a water molecule has to be present close to the γ-phosphate and adopt a near-linear attack angle (see below). In the X-ray structure,[46] the positions of water molecules are not resolved. However, in our force field MD simulations, we found a single water molecule to repeatedly and reversibly adopt such a bridging position between the γ-phosphate and E159, the carboxylate side-chain of which is located very close to Pγ (Figure ). We find a probability of ca. 15% for such a hydrolysis-competent water configuration (Table S1 and Figure S2 in the Supporting Information). In other words, water molecules frequently exchange between hydrolysis-competent and hydrolysis-noncompetent configurations, but most of the time the transporter is not in a hydrolysis-competent state. We initiated our QM/MM simulations from a representative hydrolysis-competent snapshot taken from our previous force field MD simulations.[54]

Figure summarizes the results of our QM/MM MD simulations of the first hydrolysis step. Configurations that are representative for the reactant (ES), transition (TS1), and first intermediate product (IS1) state are shown in Figure A–C, and the PMF is shown in Figure D. The reaction proceeds from the ES at a Pγ–OW distance of 2.9 Å via a free energy maximum (denoted as TS1) at 2.0 Å to the first product state (IS1), in which the Pγ–OW distance is 1.7 Å; i.e., the formation of the new bond is completed. The products of this first reaction step are ADP + HPO42–. The activation free energy is 14.1 kcal/mol, and IS1 lies 8.8 kcal/mol above the ES. These findings are confirmed by additional PMF simulations with a larger QM subsystem that also includes the side-chains of K39 and H191 (Figure D, gray line) and by PMF simulations in which the Pγ–Pβ distance was used as reaction coordinate (Figures S3 and S4 in the Supporting Information).

In the ES, the attacking water molecule forms hydrogen bonds with the E159 carboxylate group and the backbone of S163. These hydrogen bonds not only polarize the water molecule and thus activate it for proton abstraction and nucleophilic attack, but also hold it in line with the γ-phosphate (OW–Pγ–Oβ angle of about 170°). A similar geometry of the hydrolysis-competent state was proposed previously on the basis of force field MD simulations of the ABC exporter Sav1866.[55] Upon the attack of the lytic water at Pγ, the Pγ–Oβ bond breaks, and a proton is transferred from the water molecule to E159, which acts as catalytic base. This proton transfer (PT) is not imposed by the chosen RC, which is the Pγ–OW distance. Transfer of the water proton to the Oϵ atom of E159 occurs spontaneously (i.e., without any bias) at Pγ–OW distances ≤2.0 Å, i.e., at the TS and beyond (Figure E).

Figure F shows that the concerted cleavage of the Pγ–Oβ phosphate bond and PT to E159 occur in a largely synchronous fashion in the Pγ–OW distance range 1.8–2.1 Å. For larger distances, the Pγ–Oβ bond is elongated but not yet completely broken. These two modes are thus tightly coupled, underlining the importance of PT for phosphate bond cleavage. In Figure G, the ES, TS1, and IS1 trajectories are projected onto the configuration space spanned by the Pγ–Oβ distance and the distance between the E159 Oϵ atom and the water oxygen. At the TS, a broad range of Pγ–Oβ distances is sampled (magenta curve), such that the corresponding distance distribution almost overlaps with those of ES and IS1, which are much more narrow. This demonstrates that, at neighboring points along the chosen RC, there is overlap of the sampled distributions along these degrees of freedom, which are not part of the chosen RC but are of crucial importance for the reaction. Furthermore, Figure G shows that the Oϵ–OW distance is tightened at TS1, indicating that the transition state is stabilized by a strengthening of the hydrogen bond between E159 and the attacking water molecule. This observation agrees with the interpretation of X-ray crystal structures of the ABC transporter MalK trapped with transition state analogs.[36]

In addition to the attack angle, also the orientation of the lytic water molecule with respect to the γ-phosphate is important. To enable P–O bond formation, the water oxygen needs to point toward Pγ. We observed the reorientation of the water dipole in our simulations (Figure ). Upon an increase in the Pγ–OW distance from 2.9 to 3.3 Å, the free energy increases by about 1.6 kcal/mol (Figure D). This increase is due to the rearrangement of hydrogen bonds of the attacking water molecule. The hydrogen bond to the backbone of S163, which is present in the hydrolysis-competent geometry at the ES, is broken and replaced by a hydrogen bond with an Oγ atom of ATP (Figure A). These rearrangements occur at Pγ–OW distances beyond 3.3 Å; see Figure B for a representative trajectory. At smaller distances, the water molecule preferentially adopts the hydrolysis-competent orientation. This was confirmed by additional control simulations in which the PMF was calculated backward from 3.3 to 2.9 Å; no hysteresis was found.

Figure 2

Reorientation of the attacking water molecule. (A) Superposition of representative structures of the hydrolysis-competent (thick) and hydrolysis-noncompetent (transparent) states. (B) At a Pγ–OW distance of 3.3 Å, the water molecule reorients and rearranges its hydrogen bond network. The hydrogen bond between the catalytic glutamate and the water molecule is maintained.

The hydrogen bond between the lytic water molecule and the carbonyl group of the S163 backbone, which is present in the hydrolysis-competent state, is broken in the noncompetent state. As a consequence, the carbonyl group of S163 contacts a Hα atom of G130 (Figure A), suggesting that this interaction, which is likely very weak, has to be broken to stabilize the attacking water molecule in the hydrolysis-competent orientation.

Deprotonation of the Catalytic Base

After the first reaction step, E159 is protonated and hydrogen-bonded to the OH group of HPO42– (Figure C). This IS1 intermediate is further stabilized by hydrogen bonds of HPO42– with S163, H191, and K39. However, the reference pKa values of the carboxylic acid and phosphoric acid in water (ca. 4.5 and 7.2, respectively) suggest deprotonated E159 and H2PO4– as likely final reaction products. Back-transfer of the proton from E159 to the γ-phosphate requires the rearrangement of hydrogen bonds, because the protonated OH group of the γ-phosphate cannot accept a second proton. We thus considered two Oγ atoms as possible proton acceptors: the O1γ atom that coordinates the Mg2+ ion and the O3γ atom that is hydrogen-bonded to H191 (Figure ); the O2γ atom is pointing away from E159 and is at a too large distance for PT. Figure shows that PT from E159 to O1γ proceeds via a second barrier of 6.5 kcal/mol (i.e., TS2 is 15.3 kcal/mol above the initial ES) and leads to another metastable intermediate state (IS2) that is 11.8 kcal/mol above the initial ES. In IS2, the formed H2PO4– anion is stabilized by hydrogen bonds with H191 and E159, and the coordinate bond between the protonated O1γ atom and Mg2+ is elongated (Figure E), priming it for dissociation (see below).

Figure 3

Proton transfer from E159 to HPO42– and formation of H2PO4– + ADP. (A) Reactant, (B) transition, and (C) product states. (D) PMF along the distance between the E159 proton and O1γ as RC. (E) Distance of the Mg2+–O1γ coordinate bond during simulations of the reactant (IS1, black line) and product (IS2, cyan line) states.

The finding that the O1γ atom that coordinates Mg2+ accepts the proton from E159 might seem counterintuitive at first. However, recent NMR experiments on a GTPase showed that the O atom that coordinates the Mg2+ ion has the largest electron density,[56] suggesting that it is the strongest base in the PT step. The alternative PT to the O3γ atom, which is hydrogen-bonded to H191, involves a higher activation free energy barrier and is thus unlikely because of a larger distance over which the proton needs to be transferred, and because of structural strain (Figure S5 in the Supporting Information).

The described PT resets the catalytic cycle in terms of the protonation states of the involved amino acid side-chains. However, the free energy of IS2 is still rather high, +11.8 kcal/mol, suggesting that further rearrangements are necessary to yield the final product. Indeed, Figure shows that rotation of H2PO4–, such that the O2γ atom is replacing the protonated O1γ atom in the Mg2+ coordination shell, has a low activation free energy barrier of 3.2 kcal/mol and yields a final PS that is favorable in free energy, only +1.8 kcal/mol above the initial ES. We assign the low activation barrier to the weakening of the Mg2+–O1γ bond due to the previous PT step. In the final product configuration, the hydrogen bond network rearranges such that H2PO4– forms hydrogen bonds with ADP, E159, and H191 (Figure C). A proton is shared between H2PO4– and ADP, hopping back and forth between the two moieties in the PS simulations (Figure E). The populations of H2PO4– + ADP and HPO42– + H-ADP are about 75% and 25%, respectively, suggesting a free energy difference of about 1 kBT and a very low barrier between these two states.

Figure 4

Final rearrangement of hydrogen bond network and Mg2+ coordination. (A) Reactant, (B) transition, and (C) product states. (D) PMF along the distance between O2γ and Mg2+ as RC. (E) Distance between the proton on H2PO4– and ADP during simulation of the product state.

The final product, H2PO4– + ADP, is +1.8 kcal/mol higher in free energy than the initial reactant state; i.e., Δ is close to zero. This suggests that (i) the reaction is reversible up to this point, and (ii) phosphate bond cleavage does not provide a power stroke, in the sense that the free energy change that is directly associated with it cannot trigger the (partial) opening of the NBD dimer. However, one should bear in mind that, in the final product state reached in our simulations, H2PO4– and ADP are still bound to the transporter and coordinate the Mg2+ ion. Unbinding of these hydrolysis products involves reorganization of the Mg2+ coordination and individual solvation of the separated molecules. This is expected to be associated with a further lowering of the free energy. Our results thus suggest that release and solvation of ADP + P are essential for the exergonic nature of the overall reaction, a hypothesis that will be tested in future work. Our previous work[54] showed that H191 acts not only as a linchpin[36] to hold the ATP molecule in place, but also as a gatekeeper that regulates water access to the nucleotide-binding pocket. After breaking of the hydrogen bond to the γ-phosphate, H191 can adopt an outward-pointing conformation that is similar to the one observed in X-ray crystal structures of nucleotide-free BtuCD.[43,44] In this conformation, the occupancy of the binding pocket with water molecules is significantly enhanced.[54] This increased hydration of the binding pocket could be one of the initial steps on the way to the release of the reaction products, as water molecules can replace hydrogen bonds of ADP and P with the protein.

Two-Water Mechanism is Unfavorable

Our results show that, in BtuCD, a single water molecule can hydrolyze ATP, thereby transiently transferring a proton to the catalytic glutamate in a direct manner, i.e., without the help of other water molecules or amino acid side-chains. This is possible because the carboxylate group of the glutamate is located very close to the γ-phosphate. This structural feature is found in all ABC transporters and discriminates them from other NTPases, where PT along longer water wires has been suggested.[16,18,19,24,57,58] To investigate this possibility in BtuCD, we analyzed our force field MD simulations for the formation of a two-water bridge between the γ-phosphate of ATP and E159. Such a conformation was indeed observed (Figure S2 and Table S2 in the Supporting Information); its probability of ca. 14% is similar to that of the reactive conformation with one bridging water molecule (15%, see above).

In the two-water bridged configuration, the E159 side-chain is oriented away from the γ-phosphate and at a larger distance from it (the mean distance between the Cδ atom of E159 and Pγ is 5.7 and 7.3 Å for the one- and two-water bridged conformations, respectively). Starting from this two-water bridged ES, we calculated the PMF for nucleophilic water attack on Pγ, using the Pγ–OW1 distance as RC and including the second water molecule in the QM subsystem. However, unlike for the one-water mechanism described above, we did not observe PT from the lytic water molecule (WAT1) to the second, assisting water (WAT2), suggesting a barrier along this degree of freedom that cannot be overcome in the short QM/MM simulations. Therefore, we calculated a two-dimensional PMF that additionally includes the intramolecular OW–HW distance of WAT1. Figure S6 (Supporting Information) shows that PT from WAT1 via WAT2 onto E159 occurs at a Pγ–OW1 distance of 1.8 Å. The free energy barrier associated with this process is 24.9 kcal/mol and thus much higher than for the one-water mechanism. The intermediate formed by this reaction, HPO42– + ADP, is formally the same as previously but lies higher in free energy (22.8 kcal/mol). We assign this difference to a different hydrogen bond network (compare Figure A to Figure S6).

The two-water mechanism involves a significantly higher free energy barrier than the one-water mechanism (24.9 versus 15.3 kcal/mol) and is thus considered unlikely. This finding can explain the reduced ATPase activity of ATP-binding cassettes in which the catalytic glutamate is replaced by an aspartate.[59−61] This is an intriguing and up to now unexplained finding, because the E-to-D mutation does not remove the carboxylate moiety that is required as a catalytic base but merely positions it slightly further away from the γ-phosphate. Heterodimeric ABC transporters, which feature an ATPase active consensus site and an inactive degenerate site, have an aspartate residue in the latter. However, the degenerate site in ABC heterodimers bears additional noncanonical residues in the ABC signature motif, which might also contribute to the strongly reduced ATPase activity.

The above results suggest that the shorter aspartate side-chain requires a two-water bridge for nucleophilic attack of the lytic water molecule and proton shuttling to the carboxylate group. We speculate that, in this case, hydrolysis has to proceed via a mechanism that involves a higher free energy barrier and thus a slower rate. We carried out additional force field MD simulations of the E159D mutant of BtuCD-F. Figure S7 and Table S3 in the Supporting Information show that, indeed, configurations with two bridging water molecules between D159 and Pγ are frequently observed. However, the orientation of the aspartate carboxylate group and the precise arrangement of the two water molecules differ from the wild-type (E159), which might alter the actual activation energy barrier. Further studies of this topic will be the subject of future work.

Comparison to Experiments

The proposed mechanism is in agreement with, and can provide atomistic explanations for, available experimental data. The reaction rate estimated from the computed free energy profiles can be compared to experimental data for wild-type BtuCD-F from Locher and co-workers,[46,62] who reported a turnover number of about 1 ATP per second, or 0.5 s–1 per NBD monomer. Under the assumption that ATP hydrolysis is the rate-determining step,[63,64] this apparent experimental rate constant can be compared to our computations. We estimate a rate constant of ca. 7 s–1 from the highest computed activation free energy barrier of 15.3 kcal/mol, as obtained for the second reaction step (Figure ), using the Eyring equation k = pwkBTh–1 exp[−Δ‡/(RT)], where pw = 15% is the probability to actually find a hydrolysis-competent water molecule close to Pγ (see above). Within the statistical uncertainty in ΔG‡ of about ±0.8 kcal/mol (Figure D), the computational rate constants are between 2 and 27 s–1. Given that the rate constant is exponentially sensitive to the barrier height, the agreement between our simulations and experiment is very encouraging, although maybe even somewhat fortuitous.

Additional strong support for our proposed mechanism with the Walker B glutamate as catalytic base comes from site-directed mutagenesis studies of various ATP-binding cassettes.[39,40,59,60,65] Furthermore, several studies report on the pH dependence of ATPase activity. For example, Zaitseva and co-workers[63] found that the ATPase activities of the wild-type ABC dimer module of HlyB and the E-to-Q mutant, in which the Walker B glutamate is replaced by a glutamine, have similar pH dependencies with an optimum activity at around pH 7. At first sight, this finding seems to speak against the involvement of the glutamate as catalytic base. However, in a different ABC module (the TAP1 dimer), an aspartate-to-glutamate mutation in the catalytic dyad does alter the pH profile and increases ATPase activity,[61,66] underlining the important role of the Walker B glutamate. Another study,[67] which reports the pH dependence of the ATPase activity of the ABC transporter Pdr5, showed that mutation of the catalytic glutamate abolished ATPase activity, whereas mutation of the switch histidine had no strong effect, in agreement with our mechanism. Interestingly, Pdr5 has an optimal ATPase activity at pH 9,[67] clearly speaking against an essential functional role of a doubly protonated switch histidine. However, the detailed mechanistic interpretation of pH-dependent ATPase essays is complicated. For example, the NBD dimerization equilibrium is also pH-dependent, and this can affect the measured ATPase activities.[66]

Conclusions

The QM/MM free energy simulations presented in this work provide key insights into the catalytic strategy of ABC transporters. The ATP hydrolysis reaction proceeds in three steps. First, a polarized lytic water molecule attacks the ATP γ-phosphate, which—accompanied by proton transfer to the catalytic glutamate—leads to cleavage of the scissile Pγ–Oβ bond and formation of ADP + HPO42–. In the second step, the proton is transferred back from the catalytic glutamate to the γ-phosphate, yielding ADP + H2PO4– and resetting the catalytic cycle in terms of the protonation states of the involved amino acids. Third, these two reaction steps are followed by hydrogen bond rearrangements, which also involve transient changes in the coordination of the Mg2+ ion. The rate constant estimated from the computed activation free energy barrier agrees with experiments.

The described mechanism explains the catalytic activity by providing atomic-level details of the structure, dynamics, and energetics underlying the ATP hydrolysis reaction. The nucleotide-binding domains are highly conserved in sequence and structure among all ABC transporters,[6] and the mechanism is thus very likely relevant not only for BtuCD (or other type II importers), but also for a broad range of ABC transporters. However, although many residues involved in nucleotide binding and/or hydrolysis are conserved even across different NTPases,[5,6,18,68,69] the general mechanistic picture seems to be somewhat more diverse. For example, kinesin, myosin, and F1-ATPase share a similar ATP hydrolysis mechanism, in which a second assisting water molecule is involved in shuttling the proton from the lytic water molecule to the catalytic base.[70] In GTPases, diverse mechanisms have been proposed,[16] including substrate-assisted catalysis mechanisms in which the substrate itself acts as a base. These mechanisms differ from the one proposed here for ATP hydrolysis in ABC transporters. One specific feature that distinguishes ABC transporters from the previous systems is the close distance of the catalytic glutamate to the ATP γ-phosphate, and longer proton relay pathways are thus not required. We conclude that some, but not all, aspects of the hydrolysis reaction mechanism described here might be relevant also for other proteins, thus contributing to our understanding of the chemomechanical energy conversion in NTP-driven molecular machines.

Methods

The QM/MM MD simulations were started from pre-equilibrated snapshots taken from our recent force field MD study,[54] and we hence refer to this previous study for the details of the simulation setup. In brief, the complete BtuCD-F complex was simulated in a fully solvated palmitoyloleoylphosphatidylcholine (POPC) lipid bilayer. The simulations were initiated from a nucleotide-bound X-ray crystal structure (PDB entry 4FI3(46)) after replacing 5′-(β,γ-imido)triphosphate (AMP-PNP) by ATP and changing glutamine at position 159 to glutamate and cysteine at position 162 to asparagine to create the wild-type protein. The Amber ff99SB-ILDN force field[71,72] was used, in combination with the Berger lipid parameters[73,74] and the TIP4P-Ew[75] water model. For the ATP molecule in the second binding site that was not treated quantum mechanically (see below), the force field parameters of Meagher and co-workers[76] were used. The protonation states of titratable groups were assigned with PROPKA[77,78] in the presence of ATP. For the X-ray crystal structure, these computations yielded a pKa of 5.9 for the H-loop histidine; using 500 snapshots from previous 500 ns force field MD simulations of the ATP-bound state[54] yielded an average pKa of 5.2. H191 was thus modeled as singly protonated at the Nϵ atom (the cytoplasmic pH of E. coli is 7.2–7.8[79]). In additional control simulations with doubly protonated H191, the imidazole ring rotated and formed a hydrogen bond between the (protonated) Nδ–H and E159 in addition to the one between Nϵ–H and the ATP γ-phosphate. Such a conformation is not found in X-ray structures.[36]

For the QM/MM simulations, we used the standard QM/MM functionality of GROMACS (version 4.5.7)[80] in combination with a script that interfaces to GAUSSIAN09.[81] While the QM/MM routines have been an integral part of the GROMACS MD program since version 3.3, the interface script is available in the Supporting Information or for download at wwwuser.gwdg.de/~ggroenh/qmmm.html. The MPW1B95 density functional[49] was used together with a D3 dispersion correction[82] and the 6-31+G(d,p) basis set. To validate this quantum chemical method, we recalculated the reaction profile for dimethylphosphate hydrolysis, as described in the benchmark study of Ribeiro and co-workers.[50] The results (Table S4 in the Supporting Information) show that MPW1B95/6-31+G(d,p) energies closely agree with CCSD(T)/CBS values, with mean unsigned errors below 2 kcal/mol. In our QM/MM simulations of the ABC transporter, the QM subsystem consisted of the triphosphate chain of ATP; the Mg2+ ion and its two coordinating water molecules; the side-chains of S40, Q80, and E159; and one (or two) attacking water molecule(s), yielding 49 (52) QM atoms in total (Figure ; an example GAUSSIAN09 input file is contained in Supporting Information). Additional control simulations were carried out with a larger QM region that also contained the H191 and K39 side-chains (66 QM atoms). In terms of the water molecules close to the γ-phosphate, the chosen QM subsystem was complete, as there were no additional water molecules present in this ATP-binding site.[54] Electrostatic embedding of the QM subsystem into the MM point-charge surrounding[83] was used, and hydrogen link atoms were used to saturate the QM subsystem at cuts through covalent bonds across the QM/MM boundary (only C–C single bonds were cut). The partial charges of the MM atoms of the CH group next to the hydrogen link atom were set to zero, and the resulting (small) charge difference was distributed over the neighboring MM atoms such that the total charge was unchanged. Prior to MD simulation, all systems were energy-minimized (60 steepest descent steps). In the QM/MM MD simulations, the equations of motion were integrated with 1 fs time steps. An SCF convergence criterion of 10–8 Hartree was applied. For the standard QM subsystem (49 atoms), a single MD step took 6 min on a 2 × 10 core Xeon E5-2640 2.4 GHz node. NpT ensembles were simulated by coupling to a temperature bath at 300 K with a velocity rescaling[84] thermostat (τ = 0.1 ps). For constant pressure, semi-isotropic coupling was applied by separately coupling the lateral (x, y) and normal (z) directions of the periodic simulation cell to a pressure bath at 1 bar using a Berendsen barostat (τ = 2.0 ps) and compressibility 4.5 × 10–5 bar–1. All nonbonded interactions were truncated at a large cutoff of 3.5 nm, and the nonbonded pair-list was updated at every time step. The bond lengths in the MM subsystem were constrained using LINCS,[85] and SETTLE[86] was used to constrain the internal degrees of freedom of the MM water molecules. No bond-length constraints were applied to the bonds in the QM subsystem.

For the PMFs, three sets of free energy simulations were carried out. In the first set (reaction step 1), we chose the distance between the oxygen atom of the lytic water molecule, OW, and the Pγ atom of ATP as reaction coordinate (RC). To test the influence of this choice, we carried out additional PMF simulations with the Pγ–Pβ distance as RC, with very similar results (Figures S3 and S4 in the Supporting Information). In the second set (reaction step 2), the distance between the proton on E159 and an O atom of the γ-phosphate was chosen as RC. In the third set (reaction step 3), the distance between an Oγ atom and the Mg2+ ion was chosen as RC. In addition, for the alternative reaction mechanism involving two water molecules, a two-dimensional PMF was calculated. In all cases, the RC was treated with a constraint, and the force acting on the constraint was recorded in the course of the MD simulation (by using the pull code of GROMACS). Subsequent structures along the RC to initialize the MD simulations were generated by incrementing the constraint distance (a spacing of 0.1 Å was used unless otherwise noted; additional points were introduced where necessary), followed by 60 steps steepest descent energy minimization. Then, 4 ps of QM/MM MD sampling (2.5 ps for the larger QM subsystem) at 300 K was carried out at each value of the RC. Finally, the PMF can be obtained by integrating over the mean force, PMF = ∫⟨fc⟩ dr – kBT ln ⟨z(r)−1/2⟩, where r is the variable along the path, ⟨fc⟩ is the ensemble average of the constraint force at a particular point of the RC, and z(r) is the Fixman determinant of the coordinate transformation used,[87,88] which in principle has to be taken into account. However, Schlitter and co-workers[13] have shown for phosphoester hydrolysis in water that the Fixman correction due to constraining two degrees of freedom, one of which was the Pγ–OW distance, is less than 1 kcal/mol. In our case, the correction to the PMF is even smaller because (usually) only one degree of freedom was constrained. As this correction is smaller than the statistical error, it was neglected. For the two-dimensional PMF, the forces acting on the two constrained coordinates were recorded, and the PMF in the subspace of these two coordinates was obtained as described previously.[13] The first 2 ps (first 1.25 ps for the larger QM subsystem) at each distance were discarded from the analysis of fc (Figure S8 in the Supporting Information). Statistical errors were estimated using a block averaging procedure,[89] as implemented in the gmx analyze tool of GROMACS.