The multidrug transporter AcrB transports a broad range of drugs out of the cell by means of the proton-motive force. The asymmetric crystal structure of trimeric AcrB suggests a functionally rotating mechanism for drug transport. Despite various supportive forms of evidence from biochemical and simulation studies for this mechanism, the link between the functional rotation and proton translocation across the membrane remains elusive. Here, calculating the minimum free energy pathway of the functional rotation for the complete AcrB trimer, we describe the structural and energetic basis behind the coupling between the functional rotation and the proton translocation at atomic resolution. Free energy calculations show that protonation of Asp408 in the transmembrane portion of the drug-bound protomer drives the functional rotation. The conformational pathway identifies vertical shear motions among several transmembrane helices, which regulate alternate access of water in the transmembrane as well as peristaltic motions that pump drugs in the periplasm.https://doi.org/10.7554/eLife.31715.001
Bacterial multidrug resistance (MDR) is an increasing threat to current antibiotic therapy (Li et al., 2015). Resistance nodulation cell division (RND) transporters are one of the main causes of MDR in Gram-negative bacteria. These transporters pump a wide spectrum of antibiotics out of the cell by means of proton- or sodium-motive forces, conferring MDR to the bacterium when overexpressed. Understanding the mechanism of the drug efflux process is invaluable for the treatment of bacterial infections and the design of more effective drugs or inhibitors.
In Escherichia coli, the AcrA-AcrB-TolC complex is largely responsible for MDR against many lipophilic antibiotics (Du et al., 2015). This tripartite assembly spans the periplasmic space between the inner and the outer membranes of the cell, transporting drugs from the cell to the medium. TolC, an outer membrane protein, forms a generic outer membrane channel in which drugs can passively move towards the medium. AcrB is an inner membrane protein, primarily responsible for specificity towards drugs and their uptake, as well as for energy transduction. AcrA acts as an adaptor that bridges TolC and AcrB (Du et al., 2014; Wang et al., 2017).
AcrB is one of the best characterized RND transporters in both experiments and simulations, making it a prototype for studying the drug efflux mechanism of the RND family (Pos, 2009). AcrB transports drugs from the inner membrane surface or the periplasm to the TolC channel using the proton-motive force. The structure of AcrB was first solved in a three-fold symmetric form (Murakami et al., 2002), and later in an asymmetric form (Murakami et al., 2006; Seeger et al., 2006; Sennhauser et al., 2007). AcrB is a homotrimer with a triangular-prism shape, and each protomer is made up of three domains (Figure 1A): the transmembrane (TM) domain and an extensive periplasmic portion comprising the porter and the funnel domains. The TM domain transfers protons across the inner membrane down the electrochemical gradient (Figure 1C). The porter domain is made of four subdomains, PN1, PN2, PC1 and PC2 (Figure 1B), and is responsible for drug transportation. The funnel domain has a funnel-like shape with an exit pore that indirectly connects to TolC via AcrA (Du et al., 2014; Wang et al., 2017).
In the asymmetric structure, each protomer of AcrB adopts a distinct conformation that assumes one of the three functional states in the drug transport cycle: Access (Murakami et al., 2006) (or Loose [Pos, 2009; Seeger et al., 2006]), Binding (Tight) and Extrusion (Open) states, or A, B and E, respectively. The A state has entrances (an open cleft between PC1 and PC2, and a groove between TM helices, TM8 and TM9) for drug uptake. The B state has a binding pocket (called the distal binding pocket) formed by PN1, PN2 and PC1 (Figure 1B). This pocket contains many hydrophobic residues, as well as several charged and polar residues, contributing to broad substrate ‘polyspecificity’ (Takatsuka et al., 2010; Vargiu and Nikaido, 2012). The E state has a drug exit pathway to the exit pore formed by an inclined central helix (Figure 1B). The tilt of the central helix is achieved by rigid body motions of a PC2/PN1 repeat (see Figure 1B). This motion also closes the PC1/PC2 cleft, leading to contraction of the binding pocket (Pos, 2009).
These findings from the asymmetric structure and other biochemical data led to the proposal of the functional rotation mechanism (Murakami et al., 2006; Seeger et al., 2006). Hereafter, following Yao et al. (2010), we use a three letter notation to represent the state of the trimer. For example, the BEA state represents a trimeric state in which the protomer I is in the B state, the protomer II in the E state, and the protomer III in the A state (where protomers are numbered in a counterclockwise order). In the functional rotation mechanism, for example starting from the BEA state, drug efflux is coupled with a conformational transition to the EAB state. Viewed from the top, this transition is a 120 degree rotation of functional states (not a physical rotation). This mechanism is supported by both biochemical studies (Seeger et al., 2008; Takatsuka and Nikaido, 2007) and molecular dynamics (MD) simulation studies (Schulz et al., 2010; Yao et al., 2010; Zuo et al., 2015).
The three functional states satisfy the conditions required in the alternate access model (Jardetzky, 1966): the A state (and possibly also the B state) corresponds to an inward-facing state, the E state represents an outward-facing state, and the B state is an occluded state. However, despite this excellent correspondence between the structures and states, the dynamical picture of functional rotation remains unclear. Specifically, how are the conformational changes (functional rotation) and the energy source (proton translocation) synchronized and coupled dynamically? In contrast to drug efflux through small transporters, where substrate efflux and energy transduction are spatially coupled in the TM domain (Drew and Boudker, 2016), drug efflux by AcrB occurs in the periplasmic space, which is spatially separate (~50 Å apart) from the proton translocation sites in the TM domain. To understand the whole process of drug efflux of AcrB, the energy transduction mechanism between the two distant sites should be clarified, which is the purpose of this study.
The TM domain of each protomer contains twelve α-helices (TM1–12, Figure 1C). Essential residues for proton translocation are in the middle of the TM helices, D407 and D408 on TM4, K940 on TM10, and T978 on TM11. AcrB mutants in which these residues are substituted by alanine are inactive (Guan and Nakae, 2001; Takatsuka and Nikaido, 2006). In the asymmetric structure, K940 shows a distinct side chain orientation in the E state compared with that in the A and B states (Figure 1D). Specifically, K940 tilts away from D408 and towards D407 and T978 in the E state, suggesting that transient protonation of D408 contributes to the TM conformational change. Recently, all-atom MD simulation by Yamane et al. (2013) showed that the protonation of D408 stabilized the E state, whereas deprotonated D408 induced the conformation change of the TM domain toward the A state. Eicher et al., (2014) showed, by X-ray crystallography of AcrB mutants combined with MD simulations, that these conformational changes in the TM domain were related to the alternate access of water, and in turn to the functional rotation.
Although these studies provide invaluable insights into the energy transduction mechanism, the difficulty in experimental observation of the protonation states and the time-scale limitation in all-atom MD simulations preclude direct characterization of the energy transduction process. The most direct characterization would be the observation of the functional rotation in various protonation states. To achieve this, we conducted extensive MD simulations with state-of-the-art supercomputing. Using an explicit solvent/lipid all-atom model (480,074 atoms), we identified conformational transition pathways in one step of the functional rotation for the complete AcrB trimers with different protonation states.
Here, the physically most probable pathway for the functional rotation was searched with the string method (Maragliano et al., 2006; Pan et al., 2008). In this method, instead of running a long MD simulation, a pathway connecting two known structures is optimized using multiple copies of a simulation system. The pathway is represented in an important subspace (called ‘collective variables’) with discretized beads (called ‘images’). The optimization is carried out by updating the coordinates of the images through multiple simulations toward the minimum free energy pathway. To date, the string method has been successfully applied for finding the conformational transition pathways of various protein systems, such as c-Src kinase (Gan et al., 2009), myosin VI (Ovchinnikov et al., 2011), ion channels (Lev et al., 2017; Zhu and Hummer, 2010), adenylate kinase (Matsunaga et al., 2012), β2-microglobulin (Stober and Abrams, 2012), V1-ATPase (Singharoy et al., 2017), calcium pump (Das et al., 2017), and membrane transporters (Moradi et al., 2015; Moradi and Tajkhorshid, 2014). Extensive samplings around the images are necessary for accurate estimates of the energetics or free energy changes of the whole system along the pathway (Matsunaga et al., 2012; Moradi and Tajkhorshid, 2014), where free energy changes were evaluated by performing umbrella samplings around the images. In addition, to evaluate the impacts of protonation and drug binding, we conducted alchemical free energy calculations (Klimovich et al., 2015; Chipot, 2014).
By conducting these calculations, we compared the conformational pathways and energetics of two simulation systems, each with different protonation state, to clarify the causal relationship between the protonation and the functional rotation. Free energy evaluation along the pathways shows that the protonation of D408 in the B state induces the functional rotation. Structural analysis along the conformational pathway reveals that vertical shear motions in specific TM helices regulate the alternate access of water in the TM domain, as well as the peristaltic motions pumping a drug in the porter domain. These findings provide a simple and unified view of energy transduction in AcrB.
We computed the most probable conformational pathway in the functional rotation from the BEA state to the EAB state with the string method. From the observation of the crystal structure, we postulated that the protonation of D408 in the B state induces the conformational change toward the E state. This hypothesis was tested with the setup of two different simulation systems (Figure 1E): System 1 has a protonated D408 (hereafter, denoted by D408p) in protomer II, which undergoes a transition from the E to the A state. System 2 has D408p in protomer I, undergoing transition from the B to the E state. If the hypothesis is correct, then system 1 should be stable at the initial BEA state, whereas the system 2 would trigger the functional rotation. A drug (minocycline) was explicitly included in both systems. It is bound in protomer I (the B state) at first, and then transported to the exit pore during the functional rotation. No additional drug for uptake was included in order to keep the drug–AcrB interaction as simple as possible. For the string method calculation, we first generated an initial pathway in a targeted MD simulation (see Materials and methods), in which the structure of AcrB was guided from the BEA state to the EBA state by means of steering forces, simultaneously with the drug being directed toward the exit pore. The collective variables (Cartesian coordinates of Cα atoms in the porter domain and TM helices, see Materials and methods and Figure 2—figure supplement 1) used in the string method were carefully chosen according to previous simulation studies (Yamane et al., 2013; Yao et al., 2010). The drug was not included in the collective variables because it is diffusive (Schulz et al., 2010). Thirty images were used to represent the pathway.
Figure 2A shows the free energy profiles along the pathways obtained by the string method, illustrating the convergence towards the minimum free energy pathway. The impact of protonation is revealed when the converged pathways of the two systems are compared (see Figure 2B). While system 1 has an energy minimum close to the initial BEA state (at image 5), the minimum of system 2 is shifted towards the final EAB state (at image 15), indicating that the protonation state in system 2 drives the functional rotation. The structure and protonation state of the final EAB state of system 2 is identical to that of the initial BEA state of system 1, when the molecule is rotated by 120 degrees and the bound/unbound drug is ignored. This relation suggests that the release of the drug is the main factor for the increase in free energy in the late stage (images 20–30) of system 2. The molecule may need a new drug molecule to progress to the next step of rotation. Previous simulation studies (Yamane et al., 2013; Yao et al., 2010) suggested that asymmetric trimeric states are unstable without a bound drug, while the symmetric AAA state becomes most stable, avoiding the rebinding of the released drug.
We quantified the contribution from the drug–protein interactions at each image by calculating the absolute drug binding free energy to AcrB (see Materials and methods). The absolute binding free energy in system 2 is plotted in Figure 2B with the black line. The binding free energy increases rapidly at around images 10–20, prompting detailed structural analyses of the implication (Figure 2C and D). The closing of the binding pocket was monitored by the distance between the side chains of F178 and F615 (Eicher et al., 2012), and the exit gate opening by the Q124 to Y758 distance (Schulz et al., 2010), as well as the position of the drug (Figure 2C). In images 5–15, the binding pocket collapses and the exit gate opens, indicating that the minimum free energy state (image 15 in system 2) is ready to release the drug. In system 2, the initial loss of the binding free energy is compensated for by the impact of the protonation on the functional rotation (as discussed in the next subsection). In images 10–20, the drug diffuses into the gate weakly interacting with the tunnel (Figure 2C). As discussed by Schulz et al. (2010) in their targeted MD simulation study, this process may be diffusion-limited. In summary, system 2 represents the protonation state that promotes drug extrusion mediated by the functional rotation.
How does the protonation at D408 induce the conformational change in the TM domain? To address this question, we first investigated the impact of the protonation on the minimum free energy state (image 5) of system 1 by alchemically transforming it towards system 2, that is, by protonating D408 of protomer I and deprotonating D408p of protomer II (see Materials and methods). Electrostatic potential energy maps before and after the transformation of the protonation state reveal that the TM domain of protomer I is destabilized by the protonation because of a strong positive electrostatic potential bias (Figure 3A and B). This bias is caused by the surrounding positively charged residues R971 and K940, whose side-chains orient towards D408p. The relaxation of the potential bias can be seen when inspecting image 15 of system 2, the minimum free energy state (Figure 3C and Figure 3—figure supplement 1). The side-chain of R971 (on TM11) reorients downwards (towards the cytoplasm), as does D408p (on TM4) due to repulsion from K940 (on TM10). These motions involve the downshift of TM4 and TM11, whereas TM10 remains at the original position. In order to quantify this relaxation, free energy differences over the alchemical transformation (from system 1 to system 2) were evaluated at images 5 and 15 (see Materials and methods), which yielded kcal/mol at image 5, corroborating electrostatic destabilization of the protonation state in system 2, and only kcal/mol at image 15, suggesting that the protonation no longer affects the energetics at this stage.
We examined another possible protonation scheme. Here, D924 and E346, located at the periplasmic side of the TM domain, were selected for the protonation sites according to the pKa calculation by Eicher et al. (2014). The electrostatic potential map after the protonation shows a local positive change around D924, which may stabilize the system through interactions with the oxygens of the phosphate groups of lipids, and which may help to recruit water molecules for proton uptake (see the next subsection). However, the impact of the protonation at D924 is spatially localized compared to that of D408, and the change in the electrostatic potential was not propagated to the transmembrane region. The protonation at E346 scarcely changed the electrostatic potential. Therefore, these protonation sites may not have a strong impact on the functional rotation.
By monitoring the centers-of-mass for the other TM helices along the conformational pathway, we found that downshifts occur not only in TM4 and TM11 but also in TM3, TM5 and TM6 (Figure 3—figure supplement 3). Based on their analysis of crystal structures, Eicher et al. (2014) recognized two rigid assemblies in the TM domain, R1 (TM1 and TM3 to TM6) and R2 (TM7 and TM9 to TM12), which change their mutual arrangement during the transitions among B, E and A states. In terms of R1/R2 rigid body motions, the downshift motions of the TM helices can be interpreted as a vertical shear motion of R1 against R2. Tight packing of the helices within R1, and rather weak interactions on the interface between R1 and R2 may facilitate such a large rearrangement. TM11 in R2 is the exception, because it moves downwards together with R1. Previous structural studies mostly focused on lateral shear motions and rocking motions (Du et al., 2015; Eicher et al., 2014) of R1 and R2, whereas our observations clarify that protonation mainly induces a vertical shear motion. Indeed, downshifts of the specific TM helices, including TM11, can be discerned in the crystal structures when the TM helices of different states are properly aligned (Figure 3—figure supplement 4).
The next question is how D408 becomes protonated. Because proton translocation across the TM domain requires the access of water molecules to supply protons, we investigated the distribution of water molecules in protomer 1 (Figure 3D and E, and Figure 3—figure supplement 5 for other protomers). At image 5 of system 1 with D408, a water wire channel is clearly observed from periplasmic D924 (a possible transient proton binding site to facilitate proton uptake on TM10 [Eicher et al., 2014]) to the space in the middle of the helices (called the conflux space [Fischer and Kandt, 2011]) surrounded by hydrophobic residues (Figure 3D). Also, a minor water pathway comes from E346 (another possible transient proton-binding site on TM2), which merges into the conflux space. At the cytoplasmic end, the side chain of R971 breaks the channel, preventing the leakage of protons (Eicher et al., 2014).
By contrast, water distribution at image 15 of system 2 with D408p is completely different due to the vertical shear motion in the TM domain (Figure 3E). There is no water channel coming from the periplasm, resulting from the collapse of the conflux space by the downshifts of TM4 and TM11; distances between residue pairs L400–T933 and L937–F982 decrease, preventing water molecules from entering the conflux space (Figure 3E). On the other hand, the side chain of R971 opens towards cytoplasmic bulk water, exposing the protonation site and providing a possible proton release pathway. Thus, R971 is an electrostatic gate for proton diffusion to the cytoplasm, analogous to the selectivity filter of aquaporins as discussed in previous studies (Eicher et al., 2014; Eriksson et al., 2013). In summary, these two states, image 5 of system 1 and image 15 of system 2, exemplify the alternate access of water that occurs concurrently with the proton translocation.
Next, we studied energy transduction from the TM domain to the porter domain. Does D408p in the B state dynamically induce the conformational change of the porter domain? Figure 4 shows the conformational pathways of systems 1 and 2 monitored by the distances between D408 and K940 in the TM domain (see Figure 1D for crystal structure), and between PC1 and PC2 subdomains in the porter domain (see Figure 1B). As mentioned in the Introduction, the closing rigid body motion of the PC2 subdomain is related to the opening of the drug exit pathway, as well as to the contraction of the binding pocket (Figure 1B and Figure 5). The dark blue lines indicate the initial pathways generated by the targeted MD simulation. As shown in previous studies, targeted MD simulation tends to yield pathways in which local (in our case, TM domain) and global (porter domain) motions are not correlated with each other (Ovchinnikov and Karplus, 2012). Remarkably, however, the string method yielded synchronized pathways for system 2, but not for system 1. This suggests the existence of a coupling between the remote sites (the porter and TM domains) in system 2.
How are the two sites coupled? To quantify the dynamic correlations of residue pairs, we employed the mutual information analysis, which has been successfully used in analyses of allosteric coupling compared with experiments (McClendon et al., 2009). Mutual information was calculated from snapshots along the minimum free energy pathways obtained with the string method. Figure 5—figure supplement 1A plots the number of residue pairs between the TM helices and the porter domain as a function of ‘correlation’ measured by mutual information. Clearly, system 2 (broken lines in the figure) exhibits larger correlations than does system 1 (solid lines). Protomer I, among the three protomers, shows the greatest correlation, indicating that the TM domain motion of protomer I that is induced by proton translocation is harnessed for the rearrangement of the porter domain. This means that the transition from the B to the E state is the energy-consuming step (Pos, 2009).
Correlated residue pairs of system 2 are visualized in Figure 5—figure supplement 1B (Figure 5—figure supplement 2 for both systems). Here, residue pairs whose mutual information is larger than a threshold are plotted. In protomer I (B→E), high correlations were observed between TM4–TM6, TM8 and TM11, and PN1 and PC2 (and also weakly with PN2). As identified above, these TM helices (except for TM8) are involved in the downshift motions induced by D408p. To see how the TM helices and the porter domain are correlated, we visualized the structures along the pathway (Figure 5—figure supplement 1C). In particular, the coil-to-helix transition of TM8, which leverages the rigid body motion in PN1/PC2 (see Figure 1B and Figure 5), is sterically hindered by the sidechain of F459 on TM5 in the B state (Figure 5 and Figure 5—figure supplement 1C). Downshift of TM5 causes the F459 sidechain to flip, allowing the coil-to-helix transition of TM8. Thus, the downshift of TM5 reduces the barrier for the extension of TM8, expediting the PN1/PC2 rigid body motion to open the exit gate for the drug.
Protomer I also exhibits weak linkage between PN2 and several TM helices. These weak correlations may arise from PN2’s rearrangements relative to PC1 due to the shrinkage of the distal binding pocket (see Figure 1B for the binding pocket location). The downshift of the TM helices decreases contacts between the TM domain and PN2, and may accommodate the rearrangement of PN2. Previous studies indicated that a slight tilting of TM2 is coupled to PN2 motion (Du et al., 2015; Eicher et al., 2014), which is also observed here as weak correlations between TM2 and PN2.
Unlike protomer I, protomer II (E→A) has only a local correlation network between TM5 and TM8 (and weakly with TM10) and PN1 and PC2. This relatively weak correlation is associated with the reverse helix-to-coil transition of TM8, which closes the drug exit gate. Also, weak correlations can be seen in TM1, TM4, and TM10, which are consistent with the scenario described by Yamane et al. (2013) for the E→A transition: twisting motions of TM4 and TM10 interfere with TM1, PN1, and PN2. Protomer III (A→B) has rather weak correlations between PN2 and the TM helices, which may be related to the rearrangement of PN2/PC1 on distal binding pocket formation and the influence of rising TM helices on contact with the PN2 domain. These local and weak correlations in protomers II and III suggest that the minor coupling either occurs indirectly via interfacial interactions between protomers or is mediated by drug binding, not being harnessed by proton translocation.
In the light of our results and those of previous studies (Eicher et al., 2014; Seeger et al., 2008; Takatsuka and Nikaido, 2007; Yamane et al., 2013; Yao et al., 2010), we propose the cycle of the functional rotation summarized in Figure 6. Coupling of the TM and the porter domains during B→E is the energy-consuming process. Protonation of D408 in the B state drives a vertical shear motion in the TM domain, which regulates the alternate access of water (the conflux space closes, R971 opens towards the cytoplasm) as well as the motion of the porter domain. In this process, the TM8 works as a transmitter. It’s extension induces the rigid body motion of the PN2/PC1 tandem, opening the drug exit gate and shrinking the distal binding pocket. In the E→A step, the porter domain of E relaxes towards A without systematic couplings with the TM domain, perhaps driven mainly by cooperativity (interfacial interactions) between protomers. Cross-linking experiments (Seeger et al., 2008; Takatsuka and Nikaido, 2007) and coarse-grained model simulations (Yao et al., 2010) suggested that steric clashes occur between the E protomers, making two or more E-state protomers in the trimer energetically unfavorable, and making the E→A transition occur spontaneously. In the A→B step, drug that is bound in A may ‘catalyze’ (Moradi et al., 2015) formation of the distal binding pocket in the PN2/PC1 tandem. The resulting conformational change increases the number of contacts with the TM domain and drives the TM helices upwards.
This scheme is consistent with the model proposed by Pos (2009), which is based on a combination of the alternate access model and the binding change mechanism for F0F1-ATP synthase (Boyer, 1993). The present study provides evidence that energy is consumed to release the drug (i.e., the B-to-E step). This energy consumption was predicted by Pos’s model in analogy to F0F1-ATP synthase, where the energy is used to release the ATP from the beta-subunit. This study also provides a simple molecular mechanism for energy transduction from the TM domain to the porter domain in terms of the vertical shear motion in the TM domain. The same type of vertical shear motions can be seen in other transporters that employ the so-called elevator alternating-access mechanism (Drew and Boudker, 2016; Slotboom, 2014).
Finally, we discuss the possibility of other protonation schemes that may impact on functional rotation. The computational work by Eicher et al. (2014) indicated that not only D408 but also D407 is protonated in the E state, suggesting a stoichiometry of two protons per one substrate. On the other hand, pKa calculations by Yamane et al. (2013) showed that only D408 should be protonated in the E state. Furthermore, their MD simulations showed that the protonation of both D407 and D408 destabilizes the structure of the TM domain. In this study, we have double-checked our results by examining the free-energy profiles along the pathway and alchemical free energies. The consistency of two independent calculations has confirmed that the impact of transient deprotonation/protonation of D408 is certainly enough to drive functional rotation and the extrusion of a drug.
The crystal structure of the asymmetric AcrB trimer (PDB entry: 4DX5 [Eicher et al., 2012]) with bound minocycline in the distal binding pocket of protomer I (B state) was used for MD simulations. The AcrB trimer was embedded in an equilibrated lipid bilayer membrane of POPE (1-palmitoyl-2-oleoyl-sn-glycero-3-phosphoethanolamine). We estimated the TM region using TMHMM (Sonnhammer et al., 1998) and overlaid it on the equilibrated POPE membrane, then removed lipid molecules within 1.5 Å of the protein. A large hole with a diameter of ~30 Å in the TM center of AcrB should be filled with phospholipids to avoid proton leakage across the membrane. Twelve POPE molecules were manually packed into the hole (six for the upper leaflet and six for the lower leaflet). Then, the simulation box was filled with water molecules. D408 of protomer II (E state) was protonated for system 1 while D408 of protomer I (B state) was protonated for system 2. Finally, sodium ions were added to make the net charge of the system neutral. The final system comprised 480,074 atoms.
Equilibration was performed with MARBLE (Ikeguchi, 2004) and mu2lib (publicly available at http://www.mu2lib.org) using CHARMM36 for proteins (Best et al., 2012) and lipids (Klauda et al., 2010) as force-field parameters. For minocycline, the parameters compatible with CHARMM force-fields, developed by Aleksandrov and Simonson (2009), were used. Electrostatic interactions were calculated with PME (Darden et al., 1993), and the Lennard-Jones interactions were treated with the force switch algorithm (Steinbach and Brooks, 1994) (over 8–10 Å). The symplectic integrator for rigid bodies (Ikeguchi, 2004) and SHAKE (Ryckaert et al., 1977) constraint was used with MARBLE and mu2lib, respectively, using a time step of 2 fs. After an initial energy minimization, the system was gradually heated to 300 K for 1 ns. Then, the system was equilibrated for 8 ns under NPT ensemble (300 K and 1 atm). Throughout the minimization and equilibration, positional harmonic restraints were imposed to the nonhydrogen atoms of AcrB. For the minimization and first 2 ns equilibration, positional restraints were also applied to the nonhydrogen atoms of minocycline and crystallographic water molecules.
To generate the initial pathways to be used in the string method, targeted MD simulations were performed with MARBLE. Starting from the BEA trimeric state, a targeted MD simulation of 10 ns was conducted towards the EAB state. A force constant per atom of 1.0 kcal/Å2 was used for the harmonic restraint on the root-mean-square displacement (RMSD) variable measured with the nonhydrogen atoms of AcrB. The target structure (EAB state) was created by rotating the initial structure by 120 degrees around the z-axis (perpendicular to the membrane plane). As shown by Schulz et al. (2010), drug efflux towards the exit pore in the targeted MD when pulling only the atoms of AcrB was not observed in 10 ns because the drug extrusion process is diffusion limited. Thus, in order to extrude minocycline toward the exit pore, minocycline was also pulled in our simulation. Starting from the bound form in the distal binding pocket, (the RMSD variable of) minocycline was pulled to 5 Å above the exit gate (Y758) for 10 ns using the same force constant as that of AcrB. The same protocol was applied to both system 1 and 2, yielding very similar trajectories. The trajectories were post-processed for the string method calculation: snapshots were interpolated by piece-wise linear fitting in the collective variable space and 30 equidistant points (called images) were defined as the initial pathway. Image 1 corresponds to the BEA state, and image 30 represents the EAB state. Then, the all-atom coordinates (and corresponding velocities and box size) closest to each image were chosen for the initial structure of the string method.
Collective variables for the string method were carefully chosen according to previous simulation studies. We chose Cartesian coordinates of Cα atoms in the porter domain (PN1, PN2, PC1 and PC2) and selected TM helices (TM4, TM5, TM6, TM8, TM10 and a loop connecting TM5 and TM6) of all protomers (see Figure 2—figure supplement 1). The chosen porter domain residues were those used in a previous coarse-grained study, which showed that they successfully capture the essence (thermodynamic and dynamic properties) of drug transportation (Yao et al., 2010). The TM helices were chosen because a previous all-atom simulation suggested that those particular TM helices are highly mobile (Yamane et al., 2013). According to recent studies, which investigated the impact of collective variable choice on pathway accuracy, it is important to choose such mobile domains (Matsunaga et al., 2016; Pan et al., 2014). We also chose the Cγ-atoms of D407 and D408 and the Cε-atom of K940 as collective variables in order to monitor local side-chain motions (as demonstrated in Figure 4). The drug was not included in the collective variables because it is diffusive. In total, the collective variables contain 1659 atoms, consisting of 4,977 Cartesian coordinates.
String method calculations were performed with NAMD (Phillips et al., 2005) combined with in-house scripts. After defining 30 images along the initial pathway from the targeted MD trajectory, atomistic structure around each atomic structure was relaxed by 35 ns equilibration imposing positional harmonic restraints on the collective variables to its image (Figure 2—figure supplement 2). After the relaxation, the mean forces type string method (Maragliano et al., 2006) was conducted by using positional harmonic restraints with a force constant of 0.1 kcal/Å2. Mean forces were evaluated every 4 ns of simulation and images were updated according to the calculated mean force, then reparameterized to make the images equidistant from each other (also a weak smoothing operation was applied). In order to eliminate external components (translations and rotations) in the mean forces, the Cα atoms of TM helices in snapshots were fitted to the reference structure (Branduardi and Faraldo-Gómez, 2013) before mean force estimation. The reference structure was created by averaging the crystal structures of the BEA, EAB (120 degrees rotation about the z-axis from BEA) and ABE states (240 degrees rotation), thus resulting in a 3-fold symmetric structure (Zhu and Hummer, 2010). The least-squares fittings of structures were performed only in the xy-plane. z-coordinates were kept the same as in the original snapshots because the free energies of membrane proteins are symmetric only in the xy-plane. The terminal images (images 1 and 30) were kept fixed during the first 36 ns (nine cycles in the string method), and then allowed to move to relax subtle frustrations in the crystal structure that might result from detergent molecules or crystal packing. In order to avoid any drifts of the terminal images along the pathway, tangential components to the pathway were eliminated from the mean force when updating terminal images (Zhu and Hummer, 2010). Convergences of the images were monitored by RMSDs of images from those of the initial pathway (Figure 2—figure supplement 3). Also, we confirmed that terminal images converged to the same regions as those sampled by brute-force simulations (Figure 2—figure supplement 4).
In order to obtain more statistics and to evaluate free energy profiles along pathways, umbrella samplings were carried out using a single umbrella window per image. Umbrella samplings, each of 30 ns length, were performed around the images of the initial pathway, then over 30 ns around the images after the 9th cycle in the string method, over 10 ns after the 14th cycle, over 15 ns after the 19th cycle, and over 100 ns after the 24th cycle (the final converged pathway). The positional harmonic restraints with a rather weak force constant of 0.01 kcal/Å2 were used to achieve sufficient phase space overlaps between adjacent systems. Trajectory data were post-processed by MBAR (Matsunaga et al., 2012; Shirts and Chodera, 2008) and weights for the restraint-free condition were obtained. Free-energy profiles along the pathways were evaluated using the progress coordinate (Branduardi et al., 2007; Matsunaga et al., 2016) (which measures a tangential component) for pathways. To focus on the samples within a reactive tube (Ren and Vanden-Eijnden, 2005), distance metric (Branduardi et al., 2007) (which measures an orthogonal component) from pathway was introduced, and samples that deviated further than a cutoff radius of 2,000 Å2 were ignored. Changing the cutoff radius to 1500 or 3000 Å2 did not alter the results qualitatively. Statistical uncertainties (standard errors) in the free-energy profiles were estimated from block averages. For analysis of the converged pathway (Figures 2B, C and 5), only samples from those images were used in order to eliminate any biases arising from the pathways of earlier cycles.
For the simulations with NAMD under NPT ensemble (300 K and 1 atm), Langevin dynamics and Nose-Hoover Langevin piston (Feller et al., 1995; Martyna et al., 1994) were used. Long-range electrostatic interactions were calculated using PME with a grid size of <1 Å. The Lennard-Jones interactions were treated with the force switch algorithm (over 8–10 Å). SHAKE (Ryckaert et al., 1977) and SETTLE (Miyamoto and Kollman, 1992) were applied with a time step of 2 fs. A multiple time-stepping integration scheme was used with a quintic polynomial splitting function to define local force components. Calculations were performed by a hybrid OpenMP/MPI scheme on the K computer of RIKEN AICS (Yokokawa et al., 2011).
In order to quantify contributions from the drug–protein interactions to the free energy profile along the pathway, we conducted alchemical free energy calculations (Klimovich et al., 2015) (double-decoupling method) with NAMD and we evaluated the absolute drug-binding energy for each image of system 2. The last snapshots of the minimum free energy pathway of system 2 were used for the initial structure. In the calculation, minocycline was gradually decoupled with the other atoms by using a coupling parameter λ with 25 intermediate λ-states (stratifications). The simulation length for each λ-state was 1 ns. During the decoupling simulations, the same positional harmonic restraints (force constant of 0.1 kcal/Å2) as were used for the string method were imposed on the collective variables. No restraints were imposed on minocycline because it was hard to define binding sites for images after the shrinkage of the distal binding site. Free energy differences over the decoupling were evaluated by the exponential averages (Zwanzig, 1955) of the differences of the potential energies between adjacent λ states. Statistical uncertainties were estimated from block averages. The same type of decoupling simulation was performed for minocycline solvated by TIP3P water molecules. The water box size comparable to the AcrB system was used to cancel out the system size effect (Hummer et al., 1996). From two free energy differences, an absolute binding energy was calculated for each image by considering a thermodynamic cycle. The error propagation rule was applied to obtain final statistical uncertainties. To reduce computation time, absolute binding energies were calculated only for images 1, 4, 5, 7, 10–23 and 28–30.
In addition, to evaluate the impact of protonation of D408 in either the B state or the E state, we conducted another type of alchemical free energy calculation. Imposing positional harmonic restraints (force constant of 0.1 kcal/Å2) on the collective variables of image 5 of system 1 (D408p in protomer I, and D408 in protomer II), the protonation states were alchemically transformed to that of system 2 (D408 in protomer I, and D408p in protomer II). Because of the restraints on the collective variables, this calculation mainly evaluates the contribution of protonation to the free energy differences. Transformations were carried out in both forward and backward directions, with 36 intermediate λ-states for each direction, resulting in a total simulation length of 58 ns. During the transformation, the total charge of the system was kept neutral. Free energy difference of the alchemical transformation was evaluated by BAR (Bennett, 1976) (implemented in the ParseFEP toolkit [Liu et al., 2012]) and statistical uncertainty was estimated analytically. The same type of alchemical transformation was also conducted starting from the collective variables of image 15 of system 2 to the protonation state of system 1.
In order to obtain more statistics before and after the alchemical transformation, we conducted MD simulations of 50 ns length for the initial and final state, imposing the same type of positional harmonic restraints as above. Then, for the investigation of electrostatic potential energy maps, snapshots taken from these simulations were analyzed by the k-space Gaussian split Ewald method (Shan et al., 2005), and instantaneous electrostatic potentials of the smooth part were interpolated and averaged over snapshots on the fixed 3D grids. Most of the analyses were performed with in-house MATLAB scripts (publicly available at https://github.com/ymatsunaga/mdtoolbox under the BSD 3-Clause License) (Matsunaga, 2018) and MDTraj (McGibbon et al., 2015 [copy archived at https://github.com/elifesciences-publications/mdtoolbox]). Molecular figures were generated with VMD (Humphrey et al., 1996).
While the minimum free-energy path obtained by the string method provides us with the physically most probable pathway in the collective variable space, it is not useful in characterizing the dynamic properties of conformational change at the atomistic level (e.g., atoms not involved in the collective variables). Here, in order to characterize the dynamic properties of atomistic fluctuations, we generated reactive trajectories from the umbrella sampling data around the minimum free-energy path. The reactive trajectories are defined as the portions of trajectories when, after leaving the initial BEA state, it enters first the final EBA state before returning to the BEA state (E and Vanden-Eijnden, 2010). For generating the reactive trajectories, we used the property that the minimum free-energy pathway is orthogonal to the isocommittor surface (Maragliano et al., 2006). It is known that the distribution of the reactive trajectories follows the equilibrium distribution restricted on the isocommittor surface (Ren and Vanden-Eijnden, 2005). Using this property, we resampled the reactive trajectories from the umbrella sampling data according to the weights calculated by the MBAR in the previous analysis. Under discretization by 29 slices orthogonal to the pathway, 100 reactive trajectories were resampled from the umbrella sampling data. Then, mutual information of reactive trajectories was computed for pairs of residues, one from the TM and the other from the porter domains by,
where Cartesian coordinates of the center-of-masses of sidechains were used for r1 and r2 because they better capture correlated motions involving semi-rigid regions (McClendon et al., 2014; Vanwart et al., 2012) than did the dihedral angles often used in this kind of analysis (McClendon et al., 2009). p(r1) and p(r1, r2) are a probability density of r1, and a joint probability density of r1 and r2, respectively. Unlike usual correlation coefficients, the mutual information contains information about both linear and nonlinear dependences. For numerical evaluation of the mutual information, Kraskov’s algorithm (Kraskov et al., 2004), based on k-nearest neighbor distances (without binning), was used because it is practical for multi-dimensional data sets such as Cartesian coordinate data.
Molecular mechanics models for tetracycline analogsJournal of Computational Chemistry 30:243–255.https://doi.org/10.1002/jcc.21040
Efficient estimation of free energy differences from Monte Carlo dataJournal of Computational Physics 22:245–268.https://doi.org/10.1016/0021-9991(76)90078-4
Optimization of the additive CHARMM all-atom protein force field targeting improved sampling of the backbone φ, ψ and side-chain χ(1) and χ(2) dihedral anglesJournal of Chemical Theory and Computation 8:3257–3273.https://doi.org/10.1021/ct300400x
The binding change mechanism for ATP synthase--some probabilities and possibilitiesBiochimica et Biophysica Acta (BBA) - Bioenergetics 1140:215–250.https://doi.org/10.1016/0005-2728(93)90063-L
String method for calculation of minimum free-energy paths in Cartesian space in freely-tumbling systemsJournal of Chemical Theory and Computation 9:4140–4154.https://doi.org/10.1021/ct400469w
Frontiers in free-energy calculations of biological systemsWiley Interdisciplinary Reviews: Computational Molecular Science 4:71–89.https://doi.org/10.1002/wcms.1157
Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systemsThe Journal of Chemical Physics 98:10089–10092.https://doi.org/10.1063/1.464397
Conformational transitions and alternating-access mechanism in the sarcoplasmic reticulum calcium pumpJournal of Molecular Biology 429:647–666.https://doi.org/10.1016/j.jmb.2017.01.007
Shared molecular mechanisms of membrane transportersAnnual Review of Biochemistry 85:543–572.https://doi.org/10.1146/annurev-biochem-060815-014520
Structure, mechanism and cooperation of bacterial multidrug transportersCurrent Opinion in Structural Biology 33:76–91.https://doi.org/10.1016/j.sbi.2015.07.015
Transition-path theory and path-finding algorithms for the study of rare eventsAnnual Review of Physical Chemistry 61:391–420.https://doi.org/10.1146/annurev.physchem.040808.090412
Constant pressure molecular dynamics simulation: The Langevin piston methodThe Journal of Chemical Physics 103:4613–4621.https://doi.org/10.1063/1.470648
Three ways in, one way out: water dynamics in the trans-membrane domains of the inner membrane translocase AcrBProteins: Structure, Function, and Bioinformatics 79:2871–2885.https://doi.org/10.1002/prot.23122
Partial rigid-body dynamics in NPT, NPAT and NPgammaT ensembles for proteins and membranesJournal of Computational Chemistry 25:529–541.https://doi.org/10.1002/jcc.10402
Update of the CHARMM all-atom additive force field for lipids: validation on six lipid typesThe Journal of Physical Chemistry B 114:7830–7843.https://doi.org/10.1021/jp101759q
Guidelines for the analysis of free energy calculationsJournal of Computer-Aided Molecular Design 29:397–411.https://doi.org/10.1007/s10822-015-9840-9
The challenge of efflux-mediated antibiotic resistance in Gram-negative bacteriaClinical Microbiology Reviews 28:337–418.https://doi.org/10.1128/CMR.00117-14
A toolkit for the analysis of free-energy perturbation calculationsJournal of Chemical Theory and Computation 8:2606–2616.https://doi.org/10.1021/ct300242f
String method in collective variables: minimum free energy paths and isocommittor surfacesThe Journal of Chemical Physics 125:24106.https://doi.org/10.1063/1.2212942
Minimum free energy path of ligand-induced transition in adenylate kinasePLoS Computational Biology 8:e1002555.https://doi.org/10.1371/journal.pcbi.1002555
Dimensionality of collective variables for describing conformational changes of a multi-domain proteinThe Journal of Physical Chemistry Letters 7:1446–1451.https://doi.org/10.1021/acs.jpclett.6b00317
mdtoolbox, version f8eed18GitHub.
Quantifying correlations between allosteric sites in thermodynamic ensemblesJournal of Chemical Theory and Computation 5:2486–2502.https://doi.org/10.1021/ct9001812
MDTraj: a modern open library for the analysis of molecular dynamics trajectoriesBiophysical Journal 109:1528–1532.https://doi.org/10.1016/j.bpj.2015.08.015
Settle: an analytical version of the shake and rattle algorithm for rigid water modelsJournal of Computational Chemistry 13:952–962.https://doi.org/10.1002/jcc.540130805
Computational recipe for efficient description of large-scale conformational changes in biomolecular systemsJournal of Chemical Theory and Computation 10:2866–2880.https://doi.org/10.1021/ct5002285
Free energy of conformational transition paths in biomolecules: the string method and its application to myosin VIThe Journal of Chemical Physics 134:085103.https://doi.org/10.1063/1.3544209
Analysis and elimination of a bias in targeted molecular dynamics simulations of conformational transitions: application to calmodulinThe Journal of Physical Chemistry B 116:8584–8603.https://doi.org/10.1021/jp212634z
Finding transition pathways using the string method with swarms of trajectoriesThe Journal of Physical Chemistry B 112:3432–3440.https://doi.org/10.1021/jp0777059
Assessing the accuracy of two enhanced sampling methods using EGFR kinase transition pathways: the influence of collective variable choiceJournal of Chemical Theory and Computation 10:2860–2865.https://doi.org/10.1021/ct500223p
Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanesJournal of Computational Physics 23:327–341.https://doi.org/10.1016/0021-9991(77)90098-5
Functional rotation of the transporter AcrB: insights into drug extrusion from simulationsPLoS Computational Biology 6:e1000806.https://doi.org/10.1371/journal.pcbi.1000806
Engineered disulfide bonds support the functional rotation mechanism of multidrug efflux pump AcrBNature Structural & Molecular Biology 15:199–205.https://doi.org/10.1038/nsmb.1379
Gaussian split Ewald: A fast Ewald mesh method for molecular simulationThe Journal of Chemical Physics 122:54101.https://doi.org/10.1063/1.1839571
Statistically optimal analysis of samples from multiple equilibrium statesThe Journal of Chemical Physics 129:124105.https://doi.org/10.1063/1.2978177
Chemomechanical coupling in hexameric protein-protein interfaces harnesses energy within V-Type ATPasesJournal of the American Chemical Society 139:293–310.https://doi.org/10.1021/jacs.6b10744
Structural and mechanistic insights into prokaryotic energy-coupling factor transportersNature Reviews Microbiology 12:79–87.https://doi.org/10.1038/nrmicro3175
New spherical-cutoff methods for long-range forces in macromolecular simulationJournal of Computational Chemistry 15:667–683.https://doi.org/10.1002/jcc.540150702
Energetics and mechanism of the normal-to-amyloidogenic isomerization of β2-microglobulin: on-the-fly string method calculationsThe Journal of Physical Chemistry B 116:9371–9375.https://doi.org/10.1021/jp304805v
Exploring residue component contributions to dynamical network models of allosteryJournal of Chemical Theory and Computation 8:2949–2961.https://doi.org/10.1021/ct300377a
The K computer: Japanese next-generation supercomputer development projectIEEE/ACM International Symposium on Low Power Electronics and Design. pp. 371–372.
High‐temperature equation of state by a perturbation method. II. polar gasesThe Journal of Chemical Physics 23:1915–1922.https://doi.org/10.1063/1.1740604
Yibing ShanReviewing Editor; DE Shaw Research, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Energetics and Conformational Pathways of Functional Rotation in the Multidrug Transporter AcrB" for consideration by eLife. Your article has been favorably evaluated by Richard Aldrich (Senior Editor) and two reviewers, one of whom, Yibing Shan (Reviewer #1), is a member of our Board of Reviewing Editors.
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
Aiming to elucidate the detailed structural mechanism of AcrB RND transporter, this study applied the string method for conformational free-energy calculations on a supercomputing platform to identify the minimum free energy pathway of the functional rotation of AcrB. The all-atom simulations follow the conformational change of the transporter including the travel of the drug molecule (minocycline) from the binding pocket through the gate and finally release of the drug molecule, showing that the protonation of Asp408 in a transmembrane helix of the Binding protomer drives the process. This work represents potentially important progress in our understanding of the molecular mechanism of the transport process in AcrB.
Given the readership of eLife, the manuscript should include a brief and high level description of the simulation techniques (the string method, the umbrella sampling, and the alchemic free energy calculation) in the main text in terms that a reader unfamiliar with molecular dynamics simulations can relate to. What are the basic ideas and premises behind the approach? The manuscript should explain the free energy reported in Figure 2, stating clearly that this is the free energy for the whole system if that's the case.
The specific procedure should also be explained in more detail as is it is difficult to follow. It would also help if the references to the string method in the last paragraph of the Introduction were made in a more structured way – with explanations, rather than the current lumping together of a large number of very diverse papers ranging from mathematical principles all the way to practical applications.
This work only attempts to distinguish two possible schemes. In realty the process of the functional rotation is probably more complicated, possibly involving protonation of not only Asp408. For instance, Asp924 is located in a quite hydrophobic local environment and its protonation is conceivable. The proton could even be passed from Asp408 to Asp924 by the single water file connecting them. Ideally some new simulations should be performed to address this. At the very least this possibility should be more carefully discussed.
The explanation for the increase in free energy after the minima (Figure 2A) invoking drug release without uptaking new drug molecule without offering any reasoning or evidence seems casual. The reviewers suggest that the apparent high free energy might be due to the fact that any change of protonation state is not accounted for in the simulations. The authors should consider deprotonation of Asp408, which should occur at the end of the extrusion phase according to the proposed scheme (Figure 6). The authors should also consider, as aforementioned, protonation of other residues (e.g. Asp924) following Asp408 deprotonation (a single water file seems to connect Asp408 and Asp924).
Why is the drug released from the BEA complex in System 1, if the protonation of the original E protomer stabilizes the whole complex, as shown in the image 5 of Figure 2A? How is the BEA configuration changed to EAB after protonation of the second protomer as suggested by Figure 1E?
The analysis concerning the dynamics coupling of the transmembrane domains and the porter and the funnel domains (Figure 5) is too preliminary and very unintuitive. Some abstract metric of share conformation may show that the coupling is present, which is hardly surprising, but biologically it is much more interesting to discuss/demonstrate what interactions are behind the coupling.https://doi.org/10.7554/eLife.31715.021
- Yasuhiro Matsunaga
- Mitsunori Ikeguchi
- Yasuhiro Matsunaga
- Yasuhiro Matsunaga
- Tohru Terada
- Kei Moritsugu
- Hiroshi Fujisaki
- Mitsunori Ikeguchi
- Akinori Kidera
- Yasuhiro Matsunaga
- Mitsunori Ikeguchi
- Tsutomu Yamane
- Mitsunori Ikeguchi
- Tsutomu Yamane
- Mitsunori Ikeguchi
- Tsutomu Yamane
- Mitsunori Ikeguchi
- Kei Moritsugu
- Mitsunori Ikeguchi
- Akinori Kidera
- Mitsunori Ikeguchi
- Yasuhiro Matsunaga
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
We acknowledge help from Tomio Kamada and Naoyuki Miyashita in running NAMD on the K computer. We also would like to thank Yuji Sugita, Chigusa Kobayashi, Jaewoon Jung, Motoshi Kamiya, Suyong Re, Yohei Koyama, Shun Sakuraba, Luca Maragliano, and Shoji Takada for beneficial discussions. Computational resources were provided by HOKUSAI in RIKEN, and K computer in RIKEN Advanced Institute for Computational Science by the HPCI System Research project (Project ID: hp120027). This research was partly supported by Research and Development of the Next-Generation Integrated Simulation of Living Matter, by RIKEN Advanced Institute for Computational Science (to YM), by KAKENHI (Grant Numbers 24770159 [to YM] and 25291036 [to MI]), by JST PRESTO (Grant No. JPMJPR1679 to YM), by the Platform Project for Supporting Drug Discovery and Life Science Research from AMED (to KM, MI and AK), by Innovative Drug Discovery Infrastructure through Functional Control of Biomolecular Systems, by Priority Issue one in Post-K Supercomputer Development from MEXT (Project IDs hp150269, hp160223 and hp170255, to TY and MI), and by the RIKEN Dynamic Structural Biology Project (to MI).
- Yibing Shan, Reviewing Editor, DE Shaw Research, United States
- Received: September 2, 2017
- Accepted: February 12, 2018
- Version of Record published: March 6, 2018 (version 1)
© 2018, Matsunaga et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.