Dynamics of Intact MexAB-OprM Efflux Pump: Focusing on the MexA-OprM Interface

Antibiotic efflux is one of the most critical mechanisms leading to bacterial multidrug resistance. Antibiotics are effluxed out of the bacterial cell by a tripartite efflux pump, a complex machinery comprised of outer membrane, periplasmic adaptor, and inner membrane protein components. Understanding the mechanism of efflux pump assembly and its dynamics could facilitate discovery of novel approaches to counteract antibiotic resistance in bacteria. We built here an intact atomistic model of the Pseudomonas aeruginosa MexAB-OprM pump in a Gram-negative membrane model that contained both inner and outer membranes separated by a periplasmic space. All-atom molecular dynamics (MD) simulations confirm that the fully assembled pump is stable in the microsecond timescale. Using a combination of all-atom and coarse-grained MD simulations and sequence covariation analysis, we characterized the interface between MexA and OprM in the context of the entire efflux pump. These analyses suggest a plausible mechanism by which OprM is activated via opening of its periplasmic aperture through a concerted interaction with MexA.

The low resolution of the cryo-EM map (16 Å) used for fitting may lead to partial deformation of regions with high uncertainty. In addition, the presence of unwanted molecules like detergents (specially localized in AcrB and TolC transmembrane regions) can mislead the orientation and further refinement of our model. Therefore, the MDFF fitting was applied here only to the entire MexA hexamer and to the periplasmic α-helical domains of the OprM trimer (including the interface between MexA and OprM), while position restraints were applied to the entire MexB trimer and to the transmembrane βbarrel domains of the OprM trimer (see Figure 1). The latter restraints were added in order to avoid any deformation of these components due to lack of density and to preserve the conformations that are specific to the crystal structures of MexB and OprM. During MDFF fitting, the interactions between particles are described using the CHARMM force field. The cryo-EM map is converted into a potential U md , in which prohibited high energy states are allowed. A second term U ss maintains the integrity of the secondary structure elements. Thus, during the MDFF procedure, the structure is biased towards the atomistic distribution found according to the cryo-EM density. A scaling factor of 0.3 kcal mol -1 was used during the fitting, which helps to avoid unwanted overfitting to the density map that can lead to distorted structures and unstable simulations. A total of a million fitting steps were necessary to observe convergence of the fitting (Figure S1A). More details about the flexible fitting methodology can be found in the original work 8 . The above initial model of MexAB-OprM efflux pump was posteriorly embedded in a fully hydrated double 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer system, as shown in Figure 2. This system contained 1642 lipids, 45 counter ions and 315814 water molecules, generating a total of 1.261 million particles. The protein was represented using the Amber ff99SB-ILDN 9 , the POPC lipid was represented using the updated version Lipid14 of Amber 10 , and water molecules were represented with the TIP3P model 11 . We used the inflateGRO algorithm 12 for embedding proteins within the double bilayer. Initially, the hydrophobic regions were detected by the voxel-based method 12 . Transmembrane regions in the OMF and RND were independently detected and used as targets to be embedded in the bilayers, as described in the original work 12 . Afterwards, two pre-equilibrated lipid patches were aligned with the hydrophobic regions and the fully assembled pump was inserted. Water molecules were then added such that these did not overlap with lipids and proteins. Before the production runs, the system was energy minimized and pre-equilibrated for 200 ns, restraining the backbone atoms of the protein and accounting for the membrane and water relaxation. Afterwards, restraints were released and the system was allowed to evolve according to Newton's equations.

All-atom system setup of OprM in POPC lipids.
In order to study the dynamical properties of OprM, we embedded the non-gated (closed state) outer domain protein alone in a pure POPC bilayer. The initial coordinates for OprM were downloaded from the protein data bank (PDB 1WP1) 6 . Transmembrane regions in OprM were scanned as described for the case of the fully assembled pump and posteriorly inserted in a preequilibrated POPC bilayer using the inflateGRO method 12 . After protein insertion, the bilayer consisted of 694 POPC lipids. The system was solvated with 122325 water molecules and 25 counter ions were added in order to neutralize the excess of charges in the system. The backbone of the protein was positionally restrained in order to avoid any undesirable distortion of the protein during the equilibration phase (200 ns). The final system contained 416 K particles.

All-atom system setup of OprM-MexA and OprM-MexB complexes.
We also probed the MexA-dependent aperture mechanism of the periplasmic tip of OprM by performing MD simulations using two different docked structures. The first configuration made use of the closed state of the OprM (PDB 1WP1) 6 and the MexA component that was obtained from the MDFF fitting. For the second configuration, the closed OprM component was placed in direct contact with the RND inner membrane transporter MexB (PDB 2V50) 4 . Configurations were obtained using the interactive protein docking and molecular superposition software Hex 13 . Briefly, Hex is an interactive molecular graphics program for the calculation of feasible docking modes of pairs of protein molecules. Each docked molecule is modeled using 3D expansions of real orthogonal spherical polar basis functions to encode both surface shape and electrostatic charge and potential distributions. The surface of proteins is represented using a two-term surface skin plus van der Waals steric density model, whereas the electrostatic model is derived from classical electrostatic theory. The different configurations obtained after the rigid docking approach are provided in Figure S1A. These configurations were placed in a triclinic box of 8.5x8.5x21 nm dimensions and fully solvated with 49882 water molecules. Ions were inserted to fully neutralize the excess charge of the system (15 Na + ions). Before the production runs, the atoms from the backbone were positionally restrained and equilibrated for 20 ns, allowing re-ordering of the surrounding water molecules. After relaxing the system, the constraints were released. All-atom system setup of enforced OprM-MexA rotation. For the enforced rotational MD simulation, we used the docked MexA-OprM configuration ( Figure S1A); thus, the structure containing the closed OprM domain and the MDFF-derived MexA domain was inserted in a box of 8.5x8.5x21 nm and fully solvated with 50000 water molecules. To avoid unwanted rotation of the transmembrane region of OprM, position restraints were applied to the backbone atoms corresponding to the beta-barrel region (R87-L99, A110-E125, I300-T308, G322-F335) embedded in the membrane. Before the production runs, the backbone atoms of the protein were kept in these positional restraints to allow the water molecules to fully relax. Later, only the membrane-embedded regions of OprM were restrained during the production run.
The enforced rotation procedure 14 applies a torque to a concerted set of atoms and has been successfully applied to several macromolecular systems 15 . The atoms of interest (rotation group) are subject to a rotational potential V. Each atom with position ! gets assigned an equilibrium position ! ( ), rotating at a constant angular rate ω about an axis v. The given initial positions of the selected atoms (rotation group) provide the starting configuration (t=0), which sets the reference positions ! ! . The potential assuring the rotation can be expressed as: The rotation matrix Ω and the pivot u of the axis yield forces towards the reference positions ! ( ), and effectively rotating ! . In our set up, we enforced the rotation of the MexA domain with respect of the OprM structure using a force constant (rot_k0) of 400 kJ mol -2 , and a rotation rate of 0.0001 degrees ps -1 either counter-clockwise or clockwise (see the Results section in the main text).
All-atom MD simulations. Molecular simulations (except the MDFF fitting) were carried out using the GROMACS molecular simulation package version 4.6.5 16 .
Simulations were performed using a 2 fs time step. The LINCS algorithm was applied to constrain all bond lengths with a relative geometric tolerance of 10 -4 17 . Non-bonded interactions were handled using a twin-range cutoff scheme. Within a short-range cutoff of 0.9 nm, the interactions were evaluated every time step based on a pair list updated every five-time steps. The intermediate-range interactions (up to a long-range cutoff radius of 1.4 nm) were evaluated simultaneously with each pair list update and were assumed constant in between. A PME approach 18 was used to account for electrostatic interactions with a grid spacing set to 0.15 nm. Constant temperature (303 K) was maintained by weak coupling of the solvent and solute separately to a velocity-rescaling scheme 19 with a relaxation time of 1.0 ps. During equilibration, the Berendsen algorithm 20 was used to couple the system pressure at 1.0 bar through an isotropic approach with relaxation time of 1.0 ps. After this steep, the pressure of the system was scaled using a Parrinello-Rahman barostat 21,22 . Membrane systems were quenched to 303 K and coupled using a semi-isotropic pressure approach within the X and Y plane. All systems were simulated as duplicates, except the larger fully assembled pump, which was run only once. Trajectories were run for 1 us and stored every 0.02 us for posterior analysis.

Coarse-grained system setup.
In order to compute the Potential of Mean Force (PMF) of translocating antibiotic molecules through the assembled MDR tunnel, the equilibrated system from atomistic simulations was converted into a coarse-grained (CG) representation. The MARTINI model 23 was used to represent the protein, lipids, and the drug. The solvent was represented using the polarizable water model 24 . The drug was placed within the vestibule region of the assembled pump, and overlapping water molecules were removed before simulation.
The MARTINI force field contains a vast number of parameters for the simulation of proteins 25 , lipids 23 , carbohydrates 26 , etc. However, the initial set of CG parameters for Rifampicin needs to be calibrated from all-atom MD simulations. In general, for a correct parameterization the MARTINI force field relies on a 4-to-1 mapping scheme and on the reproduction of partition coefficient between different solvents 23 . To obtain the CG coordinates of Rifampicin, an all-atom set-up of the molecule was obtained using the GAFF force field 27 . The molecule was equilibrated in a water box for 0.5 us and converted to a pseudo-CG trajectory using the center of mass of the most appropriate fine-grained particles 28 : Where the vector ! !" describes the position of the pseudo-CG bead, p is the number of atoms mapped to a given coarse bead, ! is the mass of the atom j, and ! is its coordinates. From the atomistic simulation the target bonded distribution functions (bonds, angles, dihedrals) were obtained and iteratively introduced into the set of CG beads.
The coarse-grained models of antibiotic molecules were further improved by calibrating to experimental octanol-water (logP ow ) partition coefficients. A closer inspection of the MARTINI bead types revealed that many chemical groups are already carefully parameterized matching experimental and calculated partitioning coefficients. Thus, bead selection was led by the most appropriate mapping scheme. The logP ow was calculated as a difference between the solvation free energy in aqueous (∆ ) and octanol (∆ ) phase is the partitioning free energy (∆∆ ) 29 : where R and T correspond to the universal gas constant and the temperature of the system, respectively. ∆ and ∆ were calculated directly by uncoupling the nonbonded interactions of the solute with the respective solvent using the thermodynamic integration approach 30 : Here is a coupling parameter that regulates the strength of the interaction of F B (fully uncoupled) and F A (fully coupled). !" denotes the potential energy function describing the total solute-solvent interaction. The average ⋯ is taken over the MD trajectory. Calculations were performed at 25 independent points. For each individual , simulations were run for 50 ns (AA) or 100 ns (CG) respectively. After parameterization of the CG model, the obtained logP was 1.75, in good agreement with the experimental value (3.7) 31 . The final set of bonded and non-bonded terms for the coarse-grained representation are provided in the Supporting Data section.
Coarse-grained MD simulations. Briefly, the non-bonded interactions were cut off at a distance r cut of 1.2 nm. To reduce the generation of unwanted noise, we used the standard shift function of GROMACS, in which both the energy and force smoothly vanish at the cutoff distance. The global dielectric constant was adjusted to =2.5. The LJ and Coulomb potentials were shifted from r = 0.0 and r = 0.9 nm to the cutoff distance (r = 1.2 nm), respectively. The time step used to integrate the equations of motion was 25 fs. Constant temperature at 303 K was maintained by weak coupling of the solvent and membranes separately to a Berendsen 20 heat bath with a relaxation time of 1.0 ps. Pressure was controlled using the Parrinello-Rahman barostat 21,22 set at 1.0 bar.
The drug translocation potential of mean force (PMF) by MDR was calculated using the umbrella sampling approach 32 . The simulation was composed of 220 independent points spaced 1 Å apart. A restraining potential of 1000 kJ mol -1 nm -2 was applied to the center of mass of the entire drug with respect of the center of mass of the inner membrane component (MexB domain). One micro-second long simulations per point was performed along the Z reaction coordinate. PMFs were reconstructed using the weighted histogram 33 approach, and convergence was assessed by averaging the trajectories from 10 independent blocks of 0.1 us each.        Components of this tripartite efflux pump are labeled on the right. The density map on the left corresponds to the recent high-resolution (6.5 Å) cryo-EM structure 34 , and shows the TolC α-hairpins adopting a closed periplasmic aperture. The density map on the right (which was used here for MDFF-based fitting of MexAB-OprM) corresponds to the earlier low-resolution (16 Å) cryo-EM structure 1 , and shows the TolC α-hairpins adopting an open periplasmic aperture such that there is a continuous channel between the TolC trimer and AcrA hexamer. Figure S10. Mapping of the ten predicted residue pairs onto different atomistic models of the MexA-OprM interface. These models are based on the cryo-EM-fitted structures for AcrAB-TolC that describe either an (A) "adaptor-bridging" 34,35 or (B) "adaptor-wrapping" 1 interface. MexA and OprM monomers are shown as blue and red cartoons, respectively. The ten residue pairs are shown as atomic spheres, with carbon, oxygen, and nitrogen atoms colored yellow, red, and blue, respectively. (C) Average inter-Cα distances for each of the ten residue pairs over all three protomers in an "adaptor-wrapping" model of the MexA-OprM interface. Compare with the corresponding plot for the "adaptor-bridging" model shown in Figure 3B. Error bars give s.d.

V265
L39 0.0 Table S1. Predictions of covarying residue pairs between MexB and OprM by GREMLIN. Gremlin scores are expressed as probabilities that predicted residue pairs are covarying and interacting, with values ≥ 0.70 giving high-confidence predictions. All the predicted residues pairs for MexB-OprM have scores of 0, and shown are five of these pairs.

Supporting Movies
Movie S1. MexA-dependent mechanism of iris-like opening for the OprM periplasmic aperture. Residues involved during the periplasmic aperture opening of OprM are labeled with different color codes (Orange -Leu429, Magenta -Tyr413, Lime -Asp433, Cyan -Arg436, and Pink -Ala392). OprM protomers are shown in cartoons with different colors (red, blue, and yellow). (Top-left) The counter-clockwise rotation features an increase in the area enclosed by the triplet of Leu429 residues (see Results in the main text). OprM is viewed from its periplasmic aperture. (Top-right) Changes in the potential of mean force (PMF) with increasing angle rotation are shown in order to visualize the force-dependent mechanism. (Bottom) The rearrangements leading to better cogwheel-like interactions between the OprM α-hairpins with those of MexA (depicted as translucent gray cartoons). OprM is viewed from its side.