Open Access Article
Nandan
Haloi‡
a,
Archit Kumar
Vasan‡
a,
Emily J.
Geddes
bc,
Arjun
Prasanna
cd,
Po-Chao
Wen
a,
William W.
Metcalf
cd,
Paul J.
Hergenrother
bc and
Emad
Tajkhorshid
*a
aTheoretical and Computational Biophysics Group, NIH Center for Macromolecular Modeling and Bioinformatics, Beckman Institute for Advanced Science and Technology, Department of Biochemistry, Center for Biophysics and Quantitative Biology, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA. E-mail: emad@illinois.edu
bDepartment of Chemistry, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
cCarl R. Woese Institute for Genomic Biology, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
dDepartment of Microbiology, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA
First published on 14th October 2021
Antibiotic resistance of Gram-negative bacteria is largely attributed to the low permeability of their outer membrane (OM). Recently, we disclosed the eNTRy rules, a key lesson of which is that the introduction of a primary amine enhances OM permeation in certain contexts. To understand the molecular basis for this finding, we perform an extensive set of molecular dynamics (MD) simulations and free energy calculations comparing the permeation of aminated and amine-free antibiotic derivatives through the most abundant OM porin of E. coli, OmpF. To improve sampling of conformationally flexible drugs in MD simulations, we developed a novel, Monte Carlo and graph theory based algorithm to probe more efficiently the rotational and translational degrees of freedom visited during the permeation of the antibiotic molecule through OmpF. The resulting pathways were then used for free-energy calculations, revealing a lower barrier against the permeation of the aminated compound, substantiating its greater OM permeability. Further analysis revealed that the amine facilitates permeation by enabling the antibiotic to align its dipole to the luminal electric field of the porin and form favorable electrostatic interactions with specific, highly-conserved charged residues. The importance of these interactions in permeation was further validated with experimental mutagenesis and whole cell accumulation assays. Overall, this study provides insights on the importance of the primary amine for antibiotic permeation into Gram-negative pathogens that could help the design of future antibiotics. We also offer a new computational approach for calculating free-energy of processes where relevant molecular conformations cannot be efficiently captured.
The eNTRy rules provide guidelines for a drug development framework to convert Gram-positive only antibiotics to broad spectrum antibiotics through the addition of a primary amine to compounds meeting the flexibility and globularity parameters.8 This approach was used to convert 6-deoxynybomycin (6DNM),10 a natural product derivative exclusively active against Gram-positive bacteria, into an antibiotic also effective against Gram-negative pathogens by appending an amino methyl group (6DNM-NH3)8 (Fig. 1A). Compared to 6DNM, 6DNM-NH3 shows significantly enhanced accumulation in E. coli and a 64-fold increase in antibacterial activity.8 More importantly, 6DNM-NH3 has antibacterial activity against a diverse panel of multidrug resistant Gram-negative pathogens such as A. baumannii, K. pneumonia and E. cloacae.8 This same approach has been used to expand the activity of additional Gram-positive only antibiotics,11–18 including Debio-1452 and Ribocil C, into Gram-negative active derivatives, Debio-1452-NH3 and Ribocil C-PA.11,12
![]() | ||
Fig. 1 Selective permeability of antibiotics across the E. coli porin, OmpF. (A) The conversion of 6DNM to 6DNM-NH3 by introducing a primary amine enhances permeation of the compound through OmpF (each monomer colored differently).8 (B, top, left) The radius profile (calculated using HOLE39) of the crystal structure of an OmpF monomer along the membrane normal (Z-axis) relative to the membrane center of mass (C. O. M). The constriction region (CR), defined based on a pore radius <4 Å, is highlighted in orange. The CR spans z-coordinates from 1 to 8 Å. (B, top, right) Side view of an OmpF monomer, highlighting internal loop L3 in orange and the Cα atoms of acidic and basic residues within the CR in red and blue, respectively. Purple surface depicts the pore radius profile of OmpF, calculated with the program HOLE.39 (B, bottom, right) Top-down view of OmpF, depicting the electric field vector (in green dotted arrow), projected onto the membrane plane (XY), generated due to the arrangement of the charged CR residues. The electric field was evaluated using the dipole of all the water molecules within the CR (in 100 ns of equilibrium MD simulation of membrane-embedded OmpF), following a previously described protocol.25 (C) Representation of the orientation of 6DNM-NH3 in spherical coordinates. To calculate the orientation angles, a drug vector is defined between two carbon atoms (gray spheres), such that the vector points towards the primary amine group. The inclination (θ) is the angle between and the membrane normal (Z-axis); ϕ is the angle between the projections of and the electric field vector onto the XY plane. The vector for 6DNM is defined using the same two carbon atoms and the angles are calculated using the same manner. | ||
The major proteins responsible for permeating 6DNM-NH3 through the OM of E. coli were identified through porin-knockout studies of two trimeric general diffusion β-barrel OM porins: outer membrane porin F (OmpF) and OmpC.8 In this study, we use OmpF as a model system to investigate how appending a primary amine to 6DNM can significantly enhance its permeation through the OM of E. coli using molecular dynamics (MD) simulations and free energy calculations.
The main determinant for the enhanced permeation of 6DNM-NH3 is most likely the favorable interaction of the antibiotic with the porin, especially within the constriction region (CR, Fig. 1B),19 a significantly narrowed location within the pore that acts as a bottleneck for antibiotic permeation.20–31 Additionally, the CR contains many charged residues including ladders of acidic and basic residues on the L3 loop (a long, internally folded loop) and on the inner barrel wall (Fig. 1B), respectively; these residues have been suggested to interact with antibiotics to influence permeation in previous MD simulation and experimental mutagenesis studies.20–30,32–35 The arrangement of these charged residues generates a transverse electric field across the CR (Fig. 1C), and alignment of the dipole of antibiotics with this electric field has been suggested to aid in crossing the CR.24,25,29,36–38
Despite previous MD simulation studies of antibiotic permeation through OmpF,20–30,32–35,37,38 a direct comparison between aminated and amine-free antibiotics and the mechanisms by which the primary amine group aids permeation are lacking. While our previous steered molecular dynamics (SMD) simulations of 6DNM-NH3 translocation through OmpF suggested that addition of the amine group enables interactions with the negatively charged residues of the porin, potentially assisting passage of the antibiotic,8 due to the non-equilibrium nature of SMD and the limited simulation timescale, the findings from these simulations were limited, not allowing for a quantitative description of the process. Therefore, to more comprehensively characterize the molecular mechanisms facilitating permeation of aminated antibiotics, we determined the free-energy landscape associated with translocation of 6DNM and 6DNM-NH3 through OmpF using MD simulations. Calculating these free energies, however, was challenging, because adequately sampling multiple conformations of the drug in the permeation processes using equilibrium or conventional enhanced sampling MD simulations is still computationally intractable. To overcome this challenge, we developed a novel, computationally efficient algorithm, combining Monte Carlo and graph theory, to more efficiently probe the high dimensional permeation pathways of the two antibiotics through OmpF. The resulting pathways were then used for one dimensional bias exchange umbrella sampling (1D-BEUS) simulations.
From 1D-BEUS simulations, we evaluated the free-energy of permeation for 6DNM-NH3 and 6DNM which reveals a lower energetic barrier and greater permeability for the aminated antibiotic through OmpF, in accord with our previous experiments.8 Further analysis revealed that the amine facilitates permeation by enabling 6DNM-NH3 to align its dipole to the luminal electric field of the porin and by forming favorable electrostatic interactions with charged residues while hopping along the pore. The importance of these key interactions in permeation of 6DNM-NH3 was further validated by experimental mutagenesis and whole-cell accumulation assays.
To address this issue we developed a step-wise approach: (1) create a dataset of discrete poses of the antibiotic by exhaustively exploring the translational and rotational DOFs within the porin using a grid-based workflow;51 (2) use this dataset to determine multiple energetically favorable permeation trajectories from extracellular to periplasmic space using our novel Monte Carlo based pathway search (MCPS) algorithm; (3) determine the most likely permeation pathway sampled in our MCPS trajectories using Dijkstra's algorithm; and (4) use this most likely pathway to determine energetics of the permeation process using BEUS.41–45 Details of these steps are discussed in the following sections.
![]() | ||
Fig. 2 A workflow to create a dataset of discrete antibiotic poses by exhaustively exploring the translational and rotational DOFs within the pore of an OM porin and to evaluate the antibiotic–protein interaction energy (IE) for each pose. (A, left) A set of poses of 6DNM-NH3 with multiple antibiotic orientations was generated by mapping the antibiotic along its longest axis from its C. O. M. onto a Fibonacci spherical lattice (FSL).58 Each of the generated antibiotic poses was further rotated along , with an interval of Ψ. (A, right) The set of 6DNM-NH3 orientations was placed at every point of a 3-dimensional rectangular grid. During this placement, poses with clashes between the antibiotic and the protein or ring piercings were removed. The same process was applied to generate 6DNM poses within OmpF. (B) IE between the protein and poses of 6DNM-NH3 (left) or 6DNM (right), projected along the Z-coordinate, inclination (θ) and azimuthal (ϕ) of the antibiotic. | ||
First, the residues E296, D312, and D127 in the OmpF crystal structure were protonated based on previous pKa calculations.53 We note that D312 was speculated in one study to be deprotonated,54 but this suggestion was only based on indirect inference from the crystal structure. Protonating these three residues in OmpF has been a standard practice in the literature for the past decade.26,27,36,38,55–57 Then, using the center of the CR residues (Lys16, Tyr40, Arg42, Arg82, Arg132, Tyr102, Tyr106, Asp113, Met114, Leu115, Pro116, Glu117, Phe118, Gly119 and Gly120) as the origin, a 16 × 16 × 34 search grid (1 Å spacing) was defined (Fig. 2A). Copies of 6DNM-NH3 were first placed at every grid point, and then each copy was rotated by mapping a vector,
(longest axis from the C. O. M of 6DNM-NH3 to the furthest atom, the oxygen of carbonyl2 [Fig. 4B]) onto a Fibonacci spherical lattice (FSL)58 of 25 points, generating multiple orientations at each grid point. During rotation, the C. O. M of 6DNM-NH3 was set to be at the center of the FSL, while the C. O. M of the oxygen of carbonyl2 was arranged on the FSL. Each of the 25 generated orientations at each grid point was then self-rotated along
, with an interval of Ψ = 45°, creating a total of 200 poses at each grid point. We reason that this approach would evenly survey possible orientations of 6DNM-NH3. This strategy resulted in a total of 1.7 million distinct drug–protein poses.
During our grid-based orientation search process, bad contacts between protein/drug occurred in several poses. These bad contacts were either drug–protein clashes or piercings of the protein into the rings of 6DNM-NH3 (Fig. 2A). Clashes were defined if within a pose, 6DNM-NH3 had a contact distance of <1.0 Å with at least 4 heavy atoms of the protein. Ring piercings were defined by first decomposing rings into triangles and nearby bonds from protein into line segments, and then applying a geometric test to determine if line segments intersected any triangles. Poses associated with clashes and ring piercings were removed, resulting in a total of about 700
000 distinct poses.
Each pose was subjected to 200 steps of energy minimization to relax to an energetically favorable state. Minimization was performed using the generalized Born implicit solvent (GBIS) module implemented in NAMD2
59,60 to account for the effect of solvation on interactions between the drug and protein without including explicit solvent. GBIS calculates molecular electrostatics in solvent by representing water as a dielectric continuum as described by the Poisson Boltzmann equation. The elimination of explicit solvent greatly accelerates modeling efficiency. During minimization, the protein side chains and the drug were allowed to move, while protein backbone was kept fixed in order to allow the drug to relax in its environment without causing a significant conformational change in the protein. We chose to maintain the global conformation of the protein, because the X-ray structure of OmpF is already the widest state of the porin, as we have shown in a seperate study.61 Therefore, the pathway of the antibiotic through this open conformation should be the most optimum. Furthermore, when we re-ran the pose generation step without fixing the protein backbone and used our pathway search protocol (described below) to find a new pathway; we noticed a significant RMSD of L3 due to minimization alone (Fig. S7†), which is indicative of unphysical perturbation to the protein structure caused by placement of the antibiotic. Also, the drug–protein interaction energy along the new pathway is less favorable than for the original one (Fig. S7†).
59,60 (Fig. 2B). We then walk through this energy landscape along the membrane normal (Z) to determine favorable pathways connecting extracellular and periplasmic spaces while also exploring orientation space.
To do this, we first divide the Z-coordinate space from the extracellular to the periplasmic space into N 1 Å bins,
(Fig. 3A and B). Then, we assign poses to these bins, according to the Z-coordinate of the C. O. M. of the drug in each pose. That is, a pose, Si, is assigned to a bin, Bn, if
, ΔZ and Z(Si) refer to the center of Bn, the half-width of Bn (=0.5 Å) from
, and the Z-coordinate of the center of mass (C. O. M.) of the antibiotic at a particular pose, respectively. Then, a random pose (S1) from a bin closest to the extracellular space, B1, is selected as the starting point of the trajectory. Now, we need to identify a new pose in the next bin to serve as the next step in the trajectory towards the periplasmic space (decreasing Z by 1 Å). Since we are interested in modeling a chemically relevant permeation trajectory, we also need to ensure that the orientation of the antibiotic, defined by the inclination (θ) and azimuthal (ϕ) angles, does not change substantially in the next step (mechanistically unlikely), since these angular spaces are slow DOFs. To calculate the orientation angles, a drug vector
is defined (Fig. 1C), such that the vector points towards the positively charged primary amine group. Then θ is calculated as the angle between
and the membrane normal (Z-axis), while ϕ is the angle between the projections of
and the electric field vector of the protein onto the XY plane (Fig. 1C).
![]() | ||
Fig. 3 Determination of the most likely antibiotic permeation pathways through OM porins using the dataset whose generation is described in Fig. 2. (A) Flowchart describing the Monte Carlo based pathway search (MCPS) algorithm. (B) A representative MCPS trajectory of 6DNM-NH3. To run MCPS, briefly, we first discretize the Z-coordinate space into N bins to which we assign the poses generated in Fig. 2A. Then, Monte Carlo (MC) moves are used to walk through the IE landscape generated in Fig. 2B to determine favorable trajectories connecting extracellular and periplasmic spaces. In each MC move, limited change in antibiotic orientation and position are allowed. (C) Generation of a connected graph using multiple MCPS trajectories. Boltzmann weighted densities of the 2000 trajectories are projected along the Z-coordinate, inclination (θ) and azimuthal (ϕ) of the antibiotic. (D) The most likely permeation pathway for 6DNM (red) and 6DNM-NH3 (blue), evaluated using Dijkstra's algorithm.62 | ||
![]() | ||
| Fig. 4 Comparison of the energetics of permeation for 6DNM-NH3 and 6DNM through OmpF. (A) Mean (solid) and standard deviation (shaded) of the free-energy for permeation of 6DNM-NH3 (blue) and 6DNM (red). Free energies were calculated using all the MCPS-BEUS windows and projected along the Z-coordinate of the antibiotic C. O. M. (B) Mean and standard deviation of the number of hydrogen bonds between any protein residue and each functional group in the drug (colored differently according to chemical structure shown in bottom, left) of 6DNM-NH3 (top) and 6DNM (bottom, right) at different Z-positions. The total number of hydrogen bonds for each drug is colored in black. The nomenclature and color code for each electronegative atom are shown in the bottom left panel. (C) Mean and standard deviation of the drug–protein and drug–water interaction energy at different Z-positions for 6DNM-NH3 (top) and 6DNM (bottom). | ||
To transition to the next step, we compile an allowable list of poses in the next bin (B2) with θ and ϕ within a range (Δθ = 18° and Δϕ = 36°, respectively) of θ(S0) and ϕ(S0). Within the compiled list, a trial pose, Strial is chosen randomly and a move is attempted. The move is accepted with probability:
is defined as a set of poses such that:
refer to the center of
in the Z, θ, and ϕ coordinates, respectively, and ΔZ, Δθ, and Δϕ represent the range from their respective centers (Fig. 3C). Once nodes are defined, edge weights between any two nodes
are determined by calculating the negative logarithm of transition counts (each transition weighted by the Boltzmann factor of the two connected poses making up the transition):
is the set of all poses within node
, δ(Si → Sj) is 1 if a transition occurs between Si and Sj (otherwise it is 0) and Ê(Si) is the standardized energy of pose Si. The negative logarithm was taken because Dijkstra's algorithm determines the pathway with the lowest total edge weight from source to sink node. Sink and source nodes are defined as two virtual nodes, V1 (source) and V2 (sink), each connected to the set of nodes with maximum (extracellular set) and minimum (periplasmic set)
, respectively. Each connection was initially assigned an arbitrary edge weight. We then used Dijkstra's algorithm to compute the pathway connecting V1 and V2 with the lowest overall edge weight, resulting in the most likely pathway from extracellular to periplasmic space with N nodes:
, whose centers are separated by 1 Å along the Z-coordinate.
It is important to note that the search step for the most likely pathway, by definition, outputs only one most optimal pathway and might miss additional, equally optimal pathways. Therefore, to check other possible pathways, we utilized a method commonly used in transition path theory,63 where the bottleneck transition weight of the top pathway is subtracted from it in the connected graph in order to generate a new graph. Then, Dijkstra's algorithm is re-run to obtain the next optimal pathway. This protocol can be repeated multiple times to obtain additional pathways (4 additional pathways here). We observe that the majority of these new pathways adopt similar configurations to our original top pathway (Fig. S10†). Only one of these pathways is structurally distinct from the original one (termed here “the distinct pathway”), but it has a significantly lower transition weight (54.6) than the original one (230.7) (Fig. S10†), indicating its significantly smaller contribution to the permeation dynamics of the antibiotic. Furthermore, this distinct pathway has a significantly less favorable interaction energy at the region of the greatest permeation barrier (Fig. S11†), indicative of the likely unfavorable free energy profile in this region. We only used the original top pathway for further energetic characterization.
To take into account the biologically relevant configuration of the system, for each window we built a trimeric, membrane-embedded protein–drug system. To do this, we first aligned the backbone atoms of the β-barrel residues of the monomeric protein of each window to those of a single monomer of the X-ray structure of trimeric OmpF (PDB ID: 3POX).52 Then, we merged the resulting coordinates of the aligned antibiotic-monomer system with the two additional monomers of the trimeric OmpF. In the generated trimer for each window, residues E296, D312, and D127 were protonated.36,54,55 The windows were then embedded in a symmetric membrane composed of 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) lipid molecules in each leaflet generated using the Membrane Builder module of CHARMM-GUI.64 We did not use an OM composition containing lipopolysaccharides (LPS) since the membrane composition is unlikely to influence the dynamics of the drug inside the pore. Each window was solvated with TIP3P water65 and buffered in 0.15 M NaCl to generate systems containing ∼140
000 atoms with dimensions of 120 × 120 × 100 Å3.
Before performing 1D-BEUS simulations, each drug-bound, membrane-embedded trimer was minimized using 10
000 steps of the steepest descent algorithm, and then the molecular system was allowed to relax at the center of each window during a 1 ns MD simulation while the drug and heavy atoms of the protein were restrained with a force constant of 1 kcal mol−1 Å−2. This was followed by 35 ns of 1D-BEUS simulations (until the convergence of the free-energy, Fig. S13†) using the distance along the membrane normal (Z-axis) between the drug's C. O. M and the C. O. M of the antibiotic-containing monomer as the collective variable. The force constant of each window was chosen such that the exchange ratio for each window was between 0.15 and 0.3.66 The resultant force constants were 2.0 kcal mol−1 Å−2 for all windows except for windows from the entrance to exit of the CR, which had force constants of 7.0 kcal mol−1 Å−2. As shown in Fig. S12,† using these force constants resulted in good window overlap for each drug. The first 20 ns of each window were discarded as equilibration, and the rest of the trajectories were used in evaluating the free-energy. A non-parametric variation of the weighted histogram analysis method (WHAM),67 proposed by Bartels68 and implemented by Moradi and Tajkhorshid45 was used to estimate the free-energy profile from the BEUS simulations.
000 atoms with dimensions of 120 × 120 × 100 Å3. Ten replicas of each drug were then placed at the pore mouth with different initial orientations, resulting in a total of 20 independent drug–protein systems.
The systems were then minimized using the steepest descent algorithm for 10
000 steps, followed by an initial equilibration of 1 ns, during which the heavy atoms of the protein and the drug were harmonically restrained using a force constant of 1.0 kcal mol−1 Å−2. Then, a harmonic spring with a force constant of 10 kcal mol−1 Å−2 was attached to the C. O. M. of each drug, to pull it along the Z-coordinate from the pore mouth on the extracellular side to periplasmic space through a monomer of OmpF at a constant velocity of 0.5 Å ns−1. To avoid structural deformation of the protein due to the strong non-equilibrium nature of SMD, protein heavy atoms were harmonically restrained using a force constant of 1.0 kcal mol−1 Å−2.
The SMD trajectory for each drug associated with the lowest non-equilibrium work was then used to generate 1 Å windows, except for the region between the extracellular entrance and the exit of the CR (Z = −3 to 12), where 0.5 Å windows were used to ensure adequate histogram overlap. To obtain reference free energies, additional windows were added by extending the lowest-work SMD trajectory of each drug towards the water in the extracellular and periplasmic spaces such that terminal windows were at least 10 Å away from any atom of the protein. In total, 94 windows were generated along the Z-coordinate for each drug–protein system. Then, two independent 1D-BEUS simulations were performed using the same protocol described above. This protocol resulted in a good exchange ratio (between 0.15 and 0.3) for each window as suggested previously66 and good window overlap (Fig. S2†). Free energy of each drug was then evaluated using the non-parametric WHAM.45,67,68
59,60 utilizing CHARMM36m69 and CHARMM36
70 force field parameters for proteins and lipids, respectively. The force field parameters for both drugs were generated using the CHARMM General Force Field (CGenFF)71–73 with the ParamChem server. Bonded and short-range nonbonded interactions were calculated every 2 fs, and periodic boundary conditions were employed in all three dimensions. The particle mesh Ewald (PME) method74 was used to calculate long-range electrostatic interactions every 4 fs with a grid density of 1 Å−3. A force-based smoothing function was employed for pairwise nonbonded interactions at a distance of 10 Å with a cutoff of 12 Å. Pairs of atoms whose interactions were evaluated were searched and updated every 20 fs. A cutoff (13.5 Å) slightly longer than the nonbonded cutoff was applied to search for the interacting atom pairs. Constant pressure was maintained at a target of 1 atm using the Nosé–Hoover Langevin piston method.75,76 Langevin dynamics maintained a constant temperature of 310 K with a damping coefficient, γ, of 0.5 ps−1 applied to all atoms. Simulation trajectories were collected every 10 ps.
To calculate D(Z), a set of independent simulations were initiated from the same seeds used for BEUS simulations, during which the Z position of the antibiotic was restrained by means of a harmonic potential with a force constant of 10 kcal mol−1 Å−2 (as suggested by Lee et al.78). Each simulation was performed for 10 ns. From these simulations, D(Z) was calculated in the framework of an overdamped harmonic oscillator,78–80
84 donors to BW25113 recipients via conjugation with selection for exconjugants on LB agar with kanamycin (50 μg ml−1). Because the plasmids are incapable of replication in E. coli BW25113, these exconjugants arise by integration into the recipient chromosome via homologous recombination, creating partial diploids of the ompF region, with one WT and one mutant allele separated by the plasmid vector. Kanamycin-sensitive recombinants that had resolved the partial diploid state were then selected by growth on a modified LB medium that lacked NaCl and that contained 5% sucrose. The presence of the desired mutation was verified by DNA sequencing of ompF PCR products. E. coli WM8633 (lacIq, rrnB3, Δlac4787, hsdR514, Δ(araBAD)567, Δ(rhaBAD)568, rph-1, ompF(D113F)) was constructed via allele exchange using plasmid pAP002. E. coli WM8636 (lacIq, rrnB3, Δlac4787, hsdR514, Δ(araBAD)567, Δ(rhaBAD)568, rph-1, ompF(D113V)) was constructed via allele exchange using plasmid pAP003. The mutagenic plasmids carrying the ompF mutations are derivatives of pHC001A.85 These plasmids were constructed by HiFi DNA assembly (New England Biolabs, Ipswich, MA) using NotI-digested pHC001A and a mutagenized ompF allele produced by overlap-extension PCR. Plasmid pAP002 was constructed using the following primers: ompF-F: CCGGGGGATCCACTAGTTCTAGAGCGGCCGCACGTAACTGGCGTGCAAAAC; ompF-R: AAAGCTGGAGCTCCACCGCGGTGGCGGCCGCAGCGGCGGTAATGTTCTCAA; ompF(D113V)F: TTCGCGGGTCTTAAATACGCTGTTGTTGGTTCTTTCGATTACGGCCGTAA; and ompF(D113V)R: TTACGGCCGTAATCGAAAGAACCAACAACAGCGTATTTAAGACCCGCGAA. Plasmid pAP003 was constructed using ompF-F; ompF-R; ompF(D113F)F: TTCGCGGGTCTTAAATACGCTTTCGTTGGTTCTTTCGATTACGGCCGTAA; and ompF(D113F)R: TTTACGGCCGTAATCGAAAGAACCAACGAAAGCGTATTTAAGACCCGCGAA. Bold sequences in primers ompF-F and ompF-R represent homology to pHC001A required for HiFi assembly, whereas those in ompF(D113V)F, ompF(D113V)R, ompF(D113F)F, and ompF(D113F)R represent the codon changed to produce the desired mutation.
:
1 AR20 (Acros; catalogue number: 174665000)/Sigma High Temperature (Sigma–Aldrich; catalogue number: 175633), cooled to −78 °C). Bacteria were pelleted through the oil by centrifuging at 13
000 RCF for 2 min at room temperature (with the supernatant remaining above the oil). The supernatant and oil were then removed by pipetting. To lyse the samples, each pellet was resuspended in 200 μl water, and then subjected to three freeze–thaw cycle of 3 min in liquid nitrogen followed by 3 min in a water bath at 65 °C. The lysates were pelleted at 13
000 RCF for 2 min at room temperature and the supernatant was collected (180 μl). The debris were resuspended in 100 μl methanol and pelleted as before. The supernatants were removed and combined with the previous supernatants collected. Finally, the remaining debris were removed by centrifuging at 20
000 RCF for 10 min at room temperature.
Supernatants were analyzed with the QTRAP 5500 LC/MS/MS system (Sciex) in the Metabolomics Laboratory of the Roy J. Carver Biotechnology Center, University of Illinois at Urbana-Champaign. Software Analyst 1.6.2 was used for data acquisition and analysis. The 1200 Series HPLC System (Agilent Technologies) used included a degasser, an autosampler and a binary pump. The liquid chromatography separation was performed on an Agilent Zorbax SB-Aq column (4.6 mm × 50 mm; 5 μm) with mobile phase A (0.1% formic acid in water) and mobile phase B (0.1% formic acid in acetontrile). The flow rate was 0.3 ml min−1. The linear gradient was as follows: 0–3 min: 100% A; 10–15 min: 2% A; 16–20.5 min: 100% A. The autosampler was set at 15 °C. The injection volume was 1 μl. Mass spectra were acquired under positive electrospray ionization with a voltage of 5500 V. The source temperature was 450 °C. The curtain gas, ion source gas 1, and ion source gas 2 were 33, 65 and 60 psi, respectively. Multiple reaction monitoring was used for quantitation with external calibration.
Although the 1D-BEUS simulations for each drug produced apparently converging free energies (Fig. S5†), a higher permeation barrier for 6DNM-NH3 was observed compared to 6DNM (Fig. S6†), in disagreement with our experimental finding that 6DNM-NH3 reaches higher accumulation inside E. coli.8 This discrepancy between computation and experiment may arise because the initial pathways for BEUS, derived from SMD runs, may not have captured the most likely permeation pathway for the antibiotics. A proper assessment for the optimal initial pathway would require a more aggressive sampling of the orientation space that was not achieved in our SMD runs (Fig. S3†). Although a large number of SMD runs could theoretically sample orientation space comprehensively, performing an adequate number of SMD simulations for such a complex molecular system would be computationally prohibitive.
To enhance sampling of the orientation and translational DOFs of 6DNM and 6DNM-NH3, we developed a systematic and computationally efficient algorithm named Monte Carlo based pathway search (MCPS). MCPS determines multiple permeation trajectories through OmpF from the extracellular to periplasmic side using an energetic descriptor of the system. To run MCPS, we first generate a dataset containing hundreds of thousands of discrete drug–protein poses to exhaustively explore translational and rotational DOFs of the drug (Fig. 2A). This dataset is then used to construct a multidimensional drug–protein interaction energy (IE) landscape for the translation (Z-coordinate) and orientation (inclination angle, θ, and azimuthal angle, ϕ) (Fig. 2B). It is not trivial, however, to just use a rugged multi-dimensional potential energy landscape to determine the most likely pathway.88,89 Additionally, the IE landscape lacks definitive start and end states for antibiotic permeation. Therefore, we used our MCPS algorithm to walk through the IE landscape using Monte Carlo (MC) moves to determine favorable trajectories connecting extracellular and periplasmic spaces (Fig. 3). The starting point within a trajectory is randomly selected from the poses closest to the extracellular space.
To comprehensively sample all possible pathways in our defined space, we obtained 2000 MCPS trajectories that demonstrated exhaustive and convergent sampling of the conformational space (Fig. S8 and S9†). Then, we used these trajectories to build a connected graph to be used in Dijkstra's algorithm to determine the most favorable permeation pathway (Fig. 3C and D). A detailed description for each step is provided in the Methods section. Our search algorithm shows an excellent computational efficiency of ∼100 folds over 100 SMD runs (≈92 vs. 10 k node-hours; Table 1), realizing that even this many SMD runs might fail to sufficiently sample orientation space of the drug. Although estimating the overall wall clock time for calculating free energies would largely depend on the computational resources utilized, the present approach allows one to achieve a full free energy profile for an antibiotic in a few days on commonly used clusters.
| MCPS steps | Node-hour |
|---|---|
| Conformational search | 1.1 |
| Minimization and IE calculation | 90.5 |
| MCPS | 0.0042 |
| Most likely pathway calculation | 0.02 |
| Total | 91.6 |
For 6DNM-NH3, the most likely pathway derived from MCPS is largely different from the one derived from the lowest-work SMD run; however, for 6DNM, the pathways derived from either approach are similar (Fig. S15†). The MCPS-derived pathways were then used as initial seeds for two independent 1D-BEUS simulations for either drug to obtain free-energy profiles.
The free-energy profiles revealed that within the CR, both drugs face high permeation energetic barriers, similar to the SMD-derived 1D-BEUS; however, the overall free-energy barrier for 6DNM-NH3 is significantly less than 6DNM (by 2 kcal mol−1), (Fig. 4A). The barrier arises likely because: (i) the narrowness of the CR causes a loss in entropy, which has been shown to contribute to the free-energy barrier for other antibiotics' translocations through OmpF in previous MD simulation studies;20 and (ii) the loss of drug–water interactions at the CR leads to a large desolvation enthalpic penalty (Fig. 4C), which has been shown to contribute to the overall free-energy barrier in several biological processes such as protein folding,90 and protein-ligand binding.91 However, the lower free-energy barrier for 6DNM-NH3 can be attributed to the increased protein–drug interactions enabled by the positively charged amino group (amine1) (Fig. 4B and C), which offers more enthalpic compensation to the entropic and desolvation costs within the CR for the aminated compound.
We then used the MCPS-BEUS free-energy profiles combined with our computed diffusion coefficients along the pore axis to determine permeability of either drug through OmpF. Our calculations showed a 30-fold greater permeability for the aminated drug (1.1 × 10−3 cm s−1 for 6DNM-NH3 vs. 3.3 × 10−5 cm s−1 for 6DNM), substantiating the increased accumulation of 6DNM-NH3 into E. coli we observed experimentally.8
Notably, for 6DNM-NH3, a significant decrease (by ≈3 kcal mol−1) in the energetic barrier was observed in the free-energy profile derived using MCPS-BEUS compared to the one derived using SMD-BEUS. For 6DNM, the barrier calculated by the two approaches are rather similar, possibly because the initial pathways derived from these approaches are similar (Fig. S14† and S15†).
![]() | ||
| Fig. 5 Mechanism of antibiotic permeation through OmpF. Free-energy landscapes for permeation of 6DNM-NH3 and 6DNM (calculated using all respective MCPS-BEUS windows), projected along the Z-coordinate, inclination (θ) and azimuthal (ϕ) angle of the antibiotic. All free-energy values were calculated relative to the free energy in bulk solution. Mean and standard deviation of θ and ϕ of the antibiotic in the CR at (1) the entrance (Z = 9–12 Å), (2) top (Z = 5–8 Å), (3) bottom (Z = 1–4 Å), and 4) exit (Z = −3–0 Å) of the CR, along the most likely pathway are highlighted in red. | ||
![]() | ||
| Fig. 6 Key interactions along the antibiotic permeation pathway. (Top) Snapshots of both drugs at different locations along the CR, highlighting key interacting residues. (Bottom) Hydrogen bonding between the different functional groups of 6DNM-NH3 and 6DNM and protein elements at different locations along the CR. The interactions shown have hydrogen bond probability greater than 20%. Nomenclature for the antibiotic functional groups is provided in Fig. 4. Superscripts S and M represent side- and main-chain of protein residues, respectively. | ||
We find that upon entering the CR, 6DNM-NH3 orients such that its vector (
, approximately aligned with the drug's dipole moment with a magnitude of 19 Debye), aligns with the transverse electric field in the CR (θ = 94.0 ± 9.2° and ϕ = −5.1 ± 7.7°) and maintains this orientation throughout the entire CR (from entrance to exit, Fig. 5). The observed alignment of the dipole of 6DNM-NH3 with the protein's electric field may compensate the reduced entropy of the drug within the CR, as illustrated by the appearance of a local affinity site at the CR entrance (Fig. 5). This affinity site helps to stabilize 6DNM-NH3 within the CR, which could increase the probability of the drug passing through the porin. Accordingly, previous MD simulation studies have suggested that electric field alignment of the dipole of other antibiotics and the resulting affinity sites to be beneficial in permeation through OmpF.20,24,25,27 In contrast, for 6DNM, no such local affinity site was observed at the CR entrance since the dipole moment of 6DNM is significantly weaker (8 Debye) than that of 6DNM-NH3 (Fig. 5). As a result, throughout the CR, the free energy for 6DNM in both projections is calculated to be significantly greater than that of 6DNM-NH3 (Fig. 5).
The particular orientation of 6DNM-NH3 within the CR also poses amine1 to interact with acidic residues of L3 and carbonyl2 to interact with residues of the basic ladder located at the barrel wall (Fig. 5 and 6). At the CR entrance, amine1 interacts with the side chain of D121 and side/main chain of E117, located near the extracellular terminal of L3; another factor (in addition to dipole alignment) contributing to the appearance of a local affinity site at this location (Fig. 5). Within the CR, however, amine1 breaks these interactions to interact with the side chain of D113 and main chains of L115 and G119. At this location, carbonyl2 also interacts with the side chains of K16 and R42, further stabilizing the drug in the CR. Notably, the downward movement of the antibiotic towards the periplasmic space is accompanied with a drug rotation (monitored by
), while amine1 and carbonyl2 maintain similar interactions along the pore axis; this rotation is expected to help the bulky body of the drug to advance through the narrow CR. The drug exits the CR with no further reorientation; however, interactions are broken as amine1 and carbonyl2 interact with the side chains of D107 and Q60, respectively. Overall, our analysis shows that 6DNM-NH3 permeates through the CR via a stepwise mechanism as the drug hops in order for amine1 and carbonyl2 to interact with residues along the pore. In contrast, for 6DNM, significantly fewer hydrogen bonds and a completely different orientation was observed because this antibiotic lacks a primary amine and carries a weaker dipole moment (Fig. 5 and 6).
Our molecular rationalization for the increased permeation through addition of a primary amine could be generalized to the addition of other positively charged groups. In fact, we have previously shown that the addition of other positively charged nitrogen groups, namely N-alkyl guanidiniums and pyridiniums also significantly increases the accumulation of small molecules in E. coli.92 Furthermore, our mechanism could also be generalized to other general-diffusion porins of bacteria, such as OmpC (E. coli), OmpK36 (K. pneumoniae) and OmpE36 (E. cloacae), since these porins have a similar architecture to OmpF and conserve several key residues involved in the permeation process (Fig. S16†).94,95
One such simplification is to focus the calculations on a pathway that describes the most relevant DOFs during a structural transition. Identifying this transition pathway using conventional methods can still be computationally expensive; therefore, we have developed an efficient, heuristic approach using a combination of Monte Carlo simulation and graph theory. Within this approach, we first create a dataset of discrete energy minimized molecular states using a grid-based workflow.51 This dataset could be created for any system by exploring the states along all the relevant DOFs and approximately calculating the energy associated with each state. We then use the resulting dataset to determine multiple energetically favorable trajectories along the molecular process of interest using our novel Monte Carlo based pathway search (MCPS) algorithm. Finally, we determine the most likely pathway sampled in our MCPS trajectories using Dijkstra's algorithm and use this pathway to determine the free-energy underlying the molecular process with established enhanced sampling techniques such as BEUS. It is important to note that the resulting pathways could also be used to improve the quality of other enhanced sampling methods such as metadynamics and ABF.
We applied our approach to investigate antibiotic permeation through outer membrane (OM) porins of Gram-negative bacteria; a process that involves multiple slow DOFs, most importantly translation and orientation of the drug. This application is highly biomedically relevant since these pathogens are becoming increasingly resistant to currently available antibiotics, and the development of new antibiotics targeting them has been slow.1 Specifically, we investigated the molecular basis for the conversion of Gram-positive-only to broad-spectrum antibiotics with the addition of a primary amine group as observed in our previous experiments. We found that the free-energy barrier for permeation was significantly lower for the aminated derivative corresponding to a greater permeability through OM porins. Further analysis revealed that the added amino group facilitates favorable antibiotic permeation through porins by enabling the drug molecule to align its dipole to the intrinsic electric field of the porin and by forming interactions with several charged residues. The importance of these interactions, which appear to be conserved among different species, to effective permeation across the outer membrane of Gram-negative bacteria was further validated with experimental mutagenesis and whole-cell accumulation assays. In general, the structural insights offered by this study on the mechanisms by which the primary amine enhances permeation of antibiotics could help to convert other Gram-positive-only antibiotics into broad-spectrum antibiotics. Furthermore, our novel computational method could be generalized to sample better other long-timescale molecular processes involving slow DOFs.
Footnotes |
| † Electronic supplementary information (ESI) available. See DOI: 10.1039/d1sc04445a |
| ‡ These authors contributed equally to this work. |
| This journal is © The Royal Society of Chemistry 2021 |