Structural data and early simulation studies

Computational studies have been widely used to investigate efflux machinery (Athar et al., Reference Athar, Gervasoni, Catte, Basciu, Malloci, Ruggerone and Vargiu2023). Early simulations did not focus explicitly on the proton relay network, instead utilising postulated protonation states (Seeger et al., Reference Seeger, Diederichs, Eicher, Brandstätter, Schiefner, Verrey and Pos2008). Simulations by Fischer and Kandt (Reference Fischer and Kandt2011) demonstrated alternating access for water (important for mediating proton transfer) entering the TMD from the periplasm and cytoplasm in different protomer states. Further simulations by the same group displayed flexibility in the porter domain, with opening and closing of the periplasmic cleft (Fischer and Kandt, Reference Fischer and Kandt2013). For an extensive discussion of RND simulation studies, we direct the reader to a recent review (Athar et al., Reference Athar, Gervasoni, Catte, Basciu, Malloci, Ruggerone and Vargiu2023). Here we turn our attention only to simulation studies where the protonation states of the proton relay network residues have been a primary focus.

As there are two titratable residues in the relay network (D407 and D408), it is assumed that the number of protons transferred will be one or two (Alav et al., Reference Alav, Kobylka, Kuth, Pos, Picard, Blair and Bavro2021). Based on structural data, the current consensus is that both D407 and D408 are deprotonated in the A and B protomers: K940 sits between these two residues, within salt-bridging distance of both (Figure 4a,b). This relay network arrangement has been denoted the ‘engaged’ state. In the Extrusion protomer, K940 orients away from the aspartates and towards T978/N941 (‘disengaged’ state, Figure 4c). The protonation states of D407 and D408 in the Extrusion protomer remain unclear and cannot be elucidated explicitly via state-of-the-art experimental structure-determination techniques as hydrogens are not well-resolved. Protonation may be inferred from side-chain orientations, but this can yield ambiguous or conflicting results. In the Extrusion state the distance between the G403 backbone oxygen and D407 carboxylate (<4 Å, Figure 4b) suggests D407 is protonated to participate in hydrogen bonding. Simultaneously, the D407 carboxylate is also within salt-bridging distance (~4 Å) of the K940 amine nitrogen, implying D407 is deprotonated! Experiments such as carbodiimide labelling to test protonation states (Seeger et al., Reference Seeger, von Ballmoos, Verrey and Pos2009) can be a powerful investigative tool, but the bulky reagents are restricted in the narrow water wires of the TMD. Modelling and molecular dynamics (MD) simulations have been used to address this knowledge gap. These studies generally fall into two categories, those that support: single proton transfer; or transfer of two protons per conformational cycle (Figure 4c).

Figure 4. The titratable residues of AcrB adopt different configurations in different states. (a) Close-up view of key residues of the proton relay. Each image is rotated, so viewed from the same orientation. In the Access and Binding states the relay network is considered ‘engaged’, and in the Extrusion states the network is considered ‘disengaged’. Adapted from Matsunaga et al. (Reference Matsunaga, Yamane, Terada, Moritsugu, Fujisaki, Murakami, Ikeguchi and Kidera2018). (b) Inter-residue distances for residues in the proton relay network of Gram-negative bacterial RND protein structures deposited in the PDB (raw data and relevant references available in the Supplementary material). Distances were calculated as the shortest distance between side chain oxygen/nitrogen atoms, except for G403 where the backbone carbonyl oxygen was used. (c) Conformational transition from the Binding to the Extrusion state is proposed to proceed via one protonation event (D407⁻/D408^H) or two (D407^H/D408^H). Structures shown in (a) and (c) are of PDB ID: 4DX5.

Two protons

Computational support for protonation of both D407 and D408 is based largely on structures and simulations presented by Pos and colleagues (Eicher et al., Reference Eicher, Seeger, Anselmi, Zhou, Brandstätter, Verrey, Diederichs, Faraldo-Gómez and Pos2014). Eicher et al. extended the analysis of their previous crystal structure (4DX5 (Eicher et al., Reference Eicher, Cha, Seeger, Brandstätter, El-Delik, Bohnert, Kern, Verrey, Grütter, Diederichs and Pos2012)) to focus on different conformations adopted by key residues of the AcrB TMD, including D407 and D408. They compared the different structures and calculated likely protonation states of key residues using two separate approaches.

Ten independent Metropolis Monte Carlo simulations (Brooks et al., Reference Brooks, Brooks, Mackerell, Nilsson, Petrella, Roux, Won, Archontis, Bartels, Boresch, Caflisch, Caves, Cui, Dinner, Feig, Fischer, Gao, Hodoscek, Im, Kuczera, Lazaridis, Ma, Ovchinnikov, Paci, Pastor, Post, Pu, Schaefer, Tidor, Venable, Woodcock, Wu, Yang, York and Karplus2009) were conducted to evaluate protonation states of K940, R971, and all aspartate, histidine and glutamate residues. Energies were evaluated using the Poisson-Boltzmann equation. Each simulation comprised five cycles: 100,000 steps for annealing, starting at 400 K to encourage greater exploration, followed by 500,000 steps at 300 K for equilibration. Proposed changes in protonation states were accepted or rejected based on the Metropolis criterion. By averaging the protonation state over the final 500,000 steps across all 10 simulations, a global average was obtained. The periplasmic pH was assumed to be 5.5 and the cytoplasmic pH 7.5, with a dielectric constant (εp) of 4. It was noted a pH of 7.5 across the protein led to analogous results. These simulations predicted that in the E state, D408 would be protonated 100% of the time and D407 90% of the time. The study also evaluated the likelihood of protonation for the mutants K940A, R971A, D407N and D408N. When present, both D407 and D408 were predicted to be protonated in the E states. These mutants were shown previously to be functionally inactive (Guan and Nakae, Reference Guan and Nakae2001; Takatsuka and Nikaido, Reference Takatsuka and Nikaido2006): Pos and colleagues related the change in activity of these mutants to alterations in the proton relay network.

Acknowledging the dependence of this electrostatics approach on dielectric constants, free energy perturbation (FEP) was also used to re-evaluate the protonation states of key residues. All-atom simulations of the relevant AcrB protomers were used, alchemically transforming the protonation state of specific residues to calculate the change in potential energy of the system. pK_as were calculated using known side chain analogues, yielding values of 8.2 and 10.8 in the E state for D407 and D408, respectively. This agreed with their protonation probability prediction that both D407 and D408 were protonated in the E state, and matched crystal structure side chain placements (Figure 4a). It should be noted there were large errors for D407, indicating that the calculations may not be fully converged (Table 1): while FEP can be accurate, it can also be limited by sampling issues which would bias results towards the simulation starting states (Coskun et al., Reference Coskun, Chen, Clark, Lu, Harder, Wang, Friesner and Miller2022).

Table 1. Varied pKa values were calculated for D407, D408 and K940 in the studies presented in the text

Methodology and choice of starting structure can lead to varied results. The highlighted residues are those of the key E state. Also included are values calculated for the AcrB structures 2DHH and 4DX5 using PROPKA v3.5 (Olsson et al., Reference Olsson2011; Søndergaard et al., Reference Søndergaard, Olsson, Rostkowski and Jensen2011).

Eicher et al. confirmed two protonation events would create an alternating access water wire, facilitating proton transfer to and from solution. Using unbiased all-atom MD, membrane-bound AcrB was simulated using the predicted protonation states. Access to the periplasmic and cytoplasmic solvent was found to alternate with conformational state (Fischer and Kandt, Reference Fischer and Kandt2011). Overall, the evidence presented in the study is consistent with a two-proton model.

The two-proton model has also been tested on a homologue of AcrB. A recent paper by Fairweather et al. investigated substrate binding and the effect of protonation on the RND protein MtrD (Neisseria gonorrhoeae, 4MT1; Bolla et al., Reference Bolla, Su, Do, Radhakrishnan, Kumar, Long, Chou, Delmar, Lei, Rajashankar, Shafer and Yu2014; Fairweather et al., Reference Fairweather, Gupta, Chitsaz, Booth, Brown and O’Mara2021). Studying the role of protonation, united atom (UA; Schmid et al., Reference Schmid, Eichenberger, Choutko, Riniker, Winger, Mark and van Gunsteren2011) conventional MD simulations were used, changing the protonation states of D405 and D406 (AcrB D407/D408 analogues) and the presence of progesterone substrate molecules. In the UA representation, five systems were simulated with varying protonation states and number of substrates. In the absence of substrates and protonation, MtrD adopted a symmetric A/A/A state – the putative resting state (Su et al., Reference Su, Li, Gu, Takatsuka, McDermott, Nikaido and Yu2006). Protonation of both aspartates in the absence of a substrate was sufficient to begin the opening of the exit gate (distance distributions with peaks at ~8.6 and ~ 12 Å, where the closed gate is ~5–6 Å) in the same monomer but caused instability at the relay site (K948-T985 distance distribution displays three peaks) and did not induce conformational cycling in neighbouring protomers.

Substrate binding accelerated transitions towards the asymmetric structure. Protonation in the extrusion protomer and binding of a single progesterone in the proximal binding pocket of an adjacent protomer stabilised an asymmetric state (reduced fluctuations in K948-T985 and exit site distances) but did not allow complete transition on the sampled timescales. The presence of progesterone in the binding pockets of both adjacent protomers allowed further asymmetry to develop: binding of substrates in two protomers, and protonation of the final protomer was required for the functional rotation of MtrD to be observed. Binding of multiple substrates has been shown to accelerate conformational cycling in other simulation studies (Wang et al., Reference Wang, Weng and Wang2015; Matsunaga et al., Reference Matsunaga, Yamane, Terada, Moritsugu, Fujisaki, Murakami, Ikeguchi and Kidera2018). Consistent with studies on AcrB, protonation of both aspartates induced rotation within the TMD, and the most extensive water pathways were observed in the binding protomer (Fischer and Kandt, Reference Fischer and Kandt2011; Eicher et al., Reference Eicher, Seeger, Anselmi, Zhou, Brandstätter, Verrey, Diederichs, Faraldo-Gómez and Pos2014). While this study yielded results suggesting protonation of both residues (and the presence of >1 substrate) can induce conformational cycling, the authors noted the proton relay network was never fully disengaged, despite protonation of D405 and D406. Two explanations were offered: the timescale of this transition may be beyond that sampled; or disengagement of the relay network may be driven by a different proton stoichiometry, that is, the two-proton model may not be appropriate for MtrD.

One proton

In 2013, Ikeguchi and colleagues used conventional all-atom MD simulations to investigate the possible protonation states of the titratable residues of AcrB (2DHH; (Murakami et al., Reference Murakami, Nakashima, Yamashita, Matsumoto and Yamaguchi2006; Yamane et al., Reference Yamane, Murakami and Ikeguchi2013). All possible combinations of protonation states for D407/D408 were tested in the E state protomer of this asymmetric structure. While each simulation was only 100 ns, they were able to observe differences in stability of the E-state depending on the protonation pattern. Their results suggested only D408^H D407⁻ stabilised the E-state with a disengaged relay network, as determined through inter-residue distances defined by K940, D407, D408, and T978. Protonation of only D407 led to immediate rearrangement of K940: the K940-T978 hydrogen bond was disrupted, and a salt bridge was formed between K940 and E408 which became unstable after 50 ns. D408^H D407^H was also classified as unstable as K940 became dynamic, switching interactions between D408 and T978 frequently, in addition to the salt bridge between K940 and E407 becoming unstable after 90 ns. The rearrangements and subsequent instability in both systems implies that protonation of D407 is incompatible with the E-state of this crystal structure.

They further probed these findings with two D408⁻ D407⁻ simulations from snapshots taken from the D407⁻ D408^H simulation. Deprotonation resulted in immediate re-engagement of the relay network, consistent with a transition to the A state. Principal component analysis of porter domain dynamics supported these results: D407⁻D408^H remained stable in the E state, whereas D407⁻ D408⁻ displayed a structural transition from E towards A. Together these results suggested that protonation of D408 alone stabilised the E state, and its deprotonation was sufficient to induce the transition to the A state. As stated by the authors, one weakness of their approach was the limited simulation time; the MD was unbiased, and consequently the timescale sampled was insufficient to explore full conformational transitions. We note that while all possible combinations of states were simulated, only one replicate of each system was reported.

In support of their conclusions, pK_a values for the two aspartates were calculated using PROPKA (Li et al., Reference Li, Robertson and Jensen2005) and multiconformation continuum electrostatics (MCCE; Song et al., Reference Song, Mao and Gunner2009). Both methods predicted D408 to be more readily protonated than D407, and at pH 7 only D408 would be protonated. There was a disparity in the values predicted, dependent on the method and selected parameters (Table 1). The difference in values highlights the large effect dielectric constants can have on resultant pK_as (Chan et al., Reference Chan, Wilbanks, Makhatadze and Wong2012).

Jewel et al. (Reference Jewel, Liu and Dutta2017) simulated the conformational changes of apo-AcrB utilising a hybrid coarse-grain regime. AcrB (asymmetric, 2DHH (Murakami et al., Reference Murakami, Nakashima, Yamashita, Matsumoto and Yamaguchi2006) was modelled in the UA PACE forcefield (Han et al., Reference Han, Wan, Jiang and Wu2010), and surrounding lipids and water molecules were represented in the MARTINI 2 forcefield (Marrink et al., Reference Marrink, Risselada, Yefimov, Tieleman and de Vries2007). Due to its coarse-grained nature, this regime can sample longer timescales than atomistic simulations, with the caveat of reduced resolution. Similar to the methodology of Yamane et al. (Reference Yamane, Murakami and Ikeguchi2013), this study explored the different combinations of permanent protonation states for D407/D408 in the extrusion protomer: both deprotonated, one protonated, or both protonated, for a total of four systems. When both residues were deprotonated in the E protomer, the relay network approached an engaged (Access) state: K940 moves away from T978 to form salt bridges with D407 and D408 within 100 ns. Extending the simulations to the microsecond timescale, this system also displayed closing of the exit channel (Y758–Q124 decreases from ~12 Å to ~5 Å) and opening of the periplasmic cleft, again consistent with an E → A transition.

Protonation of D408 was found to maintain the orientation of the residues observed in the extrusion state of the crystal structure (K940 was always closer to T978 than D408); the exit channel remained open; and the periplasmic cleft remained closed. This contrasts with systems where D407 was protonated: in both cases, the residues of the relay network were described as “rather dynamic”, and the exit channel opens to a greater extent than that observed in the crystal structure (>15 Å vs ~ 12 Å in the crystal structure). It was noted that this may be due to the coarse-graining scheme or the relatively large time step (5 fs), but may also indicate instability of the E-state with D407^H. Taken together, these results are consistent with the atomistic study by Yamane et al. (Reference Yamane, Murakami and Ikeguchi2013): protonation of D407 destabilises the E state; protonation of D408 stabilises the E state; deprotonation of D408 drives the Extrusion to Access transition.

Following on from this, Jewel et al. (Reference Jewel, Dinh, Liu and Dutta2020) used the same hybrid coarse-graining regime to explore conformational changes induced by changes in protonation state and the presence of an indole substrate. Four systems were simulated, except this time protonation states were chosen to induce the B → E transition (rather than E → A). Three protonation schemes for the binding monomer were simulated with indole present in the binding site of that same monomer: D407⁻ D408⁻; D407⁻ D408^H; D407^H D408⁻. The fourth system had all D407/D408 residues deprotonated, and no substrate. Both systems with protonated aspartates displayed closing of the periplasmic cleft (measured by a large increase in the T676-F563 distance) and an opening of the exit channel (Q124 and Y758 distance increase to >12 Å), suggesting a B → E transition. Interestingly, the system with bound indole but no protonation of D407/D408 also displayed opening at the exit channel, though the periplasmic cleft remained open. In the absence of substrates or protonation, the final system tended towards the A/A/A state with a closed exit channel and open periplasmic cleft: the putative resting state of RND proteins (Su et al., Reference Su, Li, Gu, Takatsuka, McDermott, Nikaido and Yu2006). These results suggest that protonation of either D407 or D408 or the binding of a substrate could induce conformational cycling in the efflux pump.

Constant pH MD (CpHMD) has also been used to investigate these residues. Yue et al. (Reference Yue, Chen, Zgurskaya and Shen2017) used hybrid-solvent continuous CpHMD to investigate the protonation states and conformational dynamics of the TMD of AcrB. Only the TMD of AcrB (4DX5 (Eicher et al., Reference Eicher, Cha, Seeger, Brandstätter, El-Delik, Bohnert, Kern, Verrey, Grütter, Diederichs and Pos2012) was simulated to reduce computational cost and ensure adequate sampling. The pH replica exchange protocol (Wallace and Shen, Reference Wallace and Shen2011) was used, with 24 replicas over the pH range 1–8.5. While conventional MD assumes fixed protonation states which must be selected a priori, CpHMD allows titratable residues to be protonated/deprotonated over the course of the simulations according to the solution pH and local environment. CpHMD also enables the calculation of the pK_a values of all titratable residues simultaneously, in contrast to traditional FEP methods which yield only one pK_a at a time while fixing the protonation states of the other residues. In this case, CpHMD was coupled with replica exchange to achieve convergence of pK_a values: by exchanging systems at varying pH values, each utilising different random walks, increased sampling of the potential energy surface of the system is achieved by accelerating the crossing of energy barriers.

Using this method, several observations were made on the pH-dependent behaviour of the titratable residues. Firstly, at pH < 4.0 the native E-state is maintained, but as pH increases a new conformational state (denoted O*/E*) is sampled: at low pH the proton relay network is disengaged (K940-T978 hydrogen bonding), but at higher pH, where both D407 and D408 are deprotonated, K940 adopts an orientation more consistent with the A-state. A large-scale motion accompanies this change: a lateral rotation within the TMD resembling the E → A transition, also observed by Eicher et al. (Reference Eicher, Seeger, Anselmi, Zhou, Brandstätter, Verrey, Diederichs, Faraldo-Gómez and Pos2014).

pK_a calculations were consistent with D407⁻/D408⁻ in the A and B states (pK_as of 2.2/3.0 in A; 2.6/2.8 in B). Due to their spatial proximity and similar pK_as, protonation of these residues is strongly coupled and D407/D408 were considered as a dyad. Titration was therefore described by a stepwise model with two macroscopic pK_as: the first/second protonation events in the B state had pK_as of 3.4/2.1. Separation of 1.3 units suggests that binding of one proton is more likely than two in the B → E transition. In the E state, the pK_as of the two protonation events increase to 5.0/3.4, with D408 having the greater pK_a (D407 3.4 vs D408 5.0) and the greater increase in pK_a (+0.8 for D407 vs + 2.2 for D408). While the authors note that their pK_a values were systematically underestimated, the relative order of the values and the difference between them were considered robust: taken together the results suggest that only D408 accepts a proton in the B → E transition.

Free energy calculation methodologies were used by Matsunaga et al. (Reference Matsunaga, Yamane, Terada, Moritsugu, Fujisaki, Murakami, Ikeguchi and Kidera2018) to investigate the pathway of AcrB functional rotation. They hypothesised that D408 protonation/deprotonation drives the conformational cycling. Two systems were generated, both starting in the BEA state (4DX5 (Eicher et al., Reference Eicher, Cha, Seeger, Brandstätter, El-Delik, Bohnert, Kern, Verrey, Grütter, Diederichs and Pos2012)) with minocycline bound in the B protomer. These systems were driven towards the EAB state using targeted all-atom MD (Klauda et al., Reference Klauda, Venable, Freites, O’Connor, Tobias, Mondragon-Ramirez, Vorobyov, MacKerell and Pastor2010; Best et al., Reference Best, Zhu, Shim, Lopes, Mittal, Feig and Mackerell2012). In the first system the D408 of protomer II (E) was protonated, predicted to stabilise the initial BEA state. In the second system, D408 of protomer I (B) was protonated, predicted to induce functional rotation to the EAB state (Figure 3a) and hence stabilise the end state. Making use of the string method (Branduardi and Faraldo-Gómez, Reference Branduardi and Faraldo-Gómez2013) and umbrella sampling (Torrie and Valleau, Reference Torrie and Valleau1977), the minimum free energy pathway between the initial and final states was calculated. It was demonstrated that the local energy minimum for system 1 lay close to the initial BEA state while for system 2 the local minimum was closer to the EAB state, concluding that protonation of D408 of the B state could drive functional rotation to the E state.

To assess how the protonation state affected the free energy, system 1 was alchemically transformed towards system 2 by protonating D408 in protomer I and deprotonating D408 in protomer II. Due to restraints on collective variables (Cartesian coordinates of Cα atoms of the porter domain, TMHs accommodating the relay network, and a flexible loop between TMHs 5 and 6; side chain atoms of D407 and D408), free energy calculations were stated to evaluate the contribution of protonation to any free energy differences. By protonating protomer I (B), the free energy of the minimum was increased by 25.9 ± 0.5 kcal/mol. Electrostatic potential maps showed that the TMD of protomer I was destabilised upon protonation due to repulsion between D408^H and proximal cationic residues (K940, R971). The free energy change on transforming the protonation state of system 2 into that of system 1 was also evaluated, but this had a reduced effect (0.4 ± 6.9 kcal/mol), suggesting protonation may no longer affect energetics at this stage in the cycle.

This study further related protonation to conformational changes in the TMD. Comparing the water distribution at the local energy minima of system 1 and system 2 showed the sheer vertical motion had created alternating access to the periplasm and cytoplasm: a possible proton release pathway consistent with previous findings (Fischer and Kandt, Reference Fischer and Kandt2011). The sheer motion and TM helix tilts were also correlated to porter domain opening/closing using mutual information analysis (McClendon et al., Reference McClendon, Friedland, Mobley, Amirkhani and Jacobson2009). This study does not discount the possibility that D407 may be protonated during conformational cycling. However, it does demonstrate that singular protonation of D408 in the B state can cause functional rotation of AcrB to a stable protonated E state, and that this functional rotation is translated to the porter domain.

Conclusions & outlook

Here we have presented a series of computational studies aiming to further elucidate details of the proton transfer mechanism in RND transporters. The literature from the last decade generally favours a one-proton model in which only D408 becomes protonated on transition from the Binding to the Extrusion state, and its deprotonation is sufficient to induce cycling back to the Access state. However, inconsistencies across these studies and conflicting experimental data remain, which cannot currently be reconciled. For example, we compared the pK_as estimated by PROPKA v3.5 for 2DHH (Murakami et al., Reference Murakami, Nakashima, Yamashita, Matsumoto and Yamaguchi2006) and 4DX5 (Eicher et al., Reference Eicher, Cha, Seeger, Brandstätter, El-Delik, Bohnert, Kern, Verrey, Grütter, Diederichs and Pos2012), finding the E state D407 in particular had differing results (6.71 vs 7.86, respectively). Due to different crystallisation conditions and ligands present, it is hard to say whether the structures are within the same native ensemble or perhaps reflect different states (Best et al., Reference Best, Lindorff-Larsen, DePristo and Vendruscolo2006). In general, the pK_as shown in Table 1 highlight the difficulty in determining protonation states, as each method carries with it its own biases (Coskun et al., Reference Coskun, Chen, Clark, Lu, Harder, Wang, Friesner and Miller2022; Wilson et al., Reference Wilson, Karttunen, de Groot and Gapsys2023). It was recently demonstrated that all commonly used high throughput pK_a estimation approaches, including PROPKA and electrostatics-based methods, are prone to inaccuracies (Wei et al., Reference Wei, Hogues and Sulea2023), so it is not trivial to tease apart these details. Incorporating well-sampled ensemble states with accurate pK_as will be essential to accurately assess protonation states. Additionally, some simulation studies have shown substrate binding to accelerate or even be required for full conformational cycling (Wang et al., Reference Wang, Weng and Wang2015; Fairweather et al., Reference Fairweather, Gupta, Chitsaz, Booth, Brown and O’Mara2021). The interdependence of protonation events and substrate recognition is yet to be fully understood. Finally, it is important to note that simulated systems will inevitably simplify true biological system complexity: studies presented here have simulated the RND transporter isolated in simple symmetric phospholipid bilayers. Of particular importance is the local environment. Not only will AcrB dynamics be affected by coupling to substrates, AcrA, and TolC in vivo, but the local proton gradient across the inner membrane will affect the likelihood of protonation for these essential residues. With the maturation of increasingly powerful experimental techniques and simulation studies able to capture greater system complexity, we anticipate further exploration of this family of proteins and their functional mechanism in the future.