Machine Learning Algorithm Identifies an Antibiotic Vocabulary for Permeating Gram-Negative Bacteria

2.1. Representation of compounds

To define a representation for each compound from which we may extract a chemical vocabulary, we begin with the two-dimensional representation of a molecule as a set of atoms and bonds connecting the atoms. Using a sliding window and considering every atom in the molecule (see Fig. S1 for an example), we identify all fragments consisting of that atom plus the atoms that lie within k bonds of it for all 1 ≤ k ≤ 10 (see Fig. 2). In total, there are 22,139 different fragments comprising the training set of 595 molecules. We represent each molecule M as a N_f = 22, 139-length vector of frequencies, $\overrightarrow{\mathcal{l}}(M)=\left[f\left(x_1,M\right),\dots ,f\left(x_N,M\right)\right]$, where every entry is the number of times a particular fragment appears (may be 0), normalized by the number of atoms in molecule M, L(M), such that,

in which n(x_i, M) represents the number of times fragment x_i appears in molecule M. Intuitively, containing larger numbers of relevant fragments should be correlated with increased activity; however, we wanted a metric that was independent of molecular size, thus our division by L(M). Although this is a somewhat naive representation akin to one-hot encoding in traditional ML³⁶, we employ it as the simplest way to test whether information about permeation ability is contained within the fragments, although greater sophistication could be achieved through the use of more complex representations^37,38.

2.2. Compound activity data

The two major contributers to the impermeability of Gram-negative bacteria in general and P. aeruginosa in specific are the OM and the efflux pumps that actively remove molecules from the periplasm and cytoplasm^2,40. To separate the effects of the efflux pumps from the effects of the OM, we have recently created different mutant strains of Gram-negative bacteria⁴¹. In this study, we focused on the effects of the OM alone by using two strategically designed mutant strains lacking the effects of efflux. In the first strain, compounds are impeded by the OM barrier, while in the second strain, they are not. Specifically, we studied mutants of the P. aeruginosa PAO1 strain. The “PΔ6” mutant is a variant of P. aeruginosa in which the genes encoding for the six best characterized efflux pumps have been deleted, which essentially removes the contribution of active efflux in antibacterial activities of antibiotics. It has no other effects; indeed, we have recently shown that there is no significant membrane disorganization introduced by deletions^8,14. The “Pore” mutant is a variant-not studied in the work-modified to contain large (~2.4 nm in diameter) pores that allow nondiscriminate entry of drugs, which essentially removes the effects of the impermeable outer membrane with no other effects on cell physiology. The “PΔ6-Pore” mutant is a variant combining both previous modifications. In this study, we focus on the PΔ6 and PΔ6-Pore mutants, which both lack efflux pumps. For the drug property input to the algorithm, we experimentally measured the MICs of over 500 compounds exhibiting antibacterial activities in at least one out of the two different mutant strains of P. aeruginosa PAO1 (see Sec. 2.2.1 for a complete description of the curated dataset). We then computed the ratio of compound MIC values in the PΔ6-Pore mutant of P. aeruginosa PAO1 to their MIC values in the PΔ6 mutant of P. aeruginosa PAO1 $\left(\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}\right)$. When this ratio goes to one, the drug’s efficacy is the same whether or not the outer membrane is intact, because the main difference between the strains is whether or not the outer membrane has been hyperporinated. We define a drug possessing this property as being a “good permeator” or saying that it “permeates well.” More specifically, we define a set of five compound classes based on the ratio $\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}$, where μ_PΔ6-Pore is the MIC of the compound in the PΔ6-Pore mutant and μ_PΔ6 is the MIC of the compound in the PΔ6 mutant. If $\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}<0.2$, $\text{Class}\left(\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}\right)=0$ (“non-permeators”); if $0.2\le \frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}<0.4$, $\text{Class}\left(\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}\right)=1$; if $0.4\le \frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}<0.6$, $\text{Class}\left(\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}\right)=2$; if $0.6\le \frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}<0.8$, $\text{Class}\left(\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}\right)=3$; and if $\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}>0.8$, $\text{Class}\left(\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}\right)=4$ (“good permeators”). The class breakdown is as follows: ≈ 48% of MIC ratios fall into class 0, 10% into class 1, 9% into class 2, 10% into class 4, and 22% into class 4.

P. aeruginosa cells were grown in Luria Bertani Broth (LB) (10 g tryptone, 5 g yeast extract, 5 g NaCl per liter, pH 7.0) at 37°C with shaking. Minimum inhibitory concentration (MIC) determination was carried out using the 2-fold broth dilution method as described previously⁸. Two independent experiments were carried out, giving precision =/-two fold. The expression of the Pore was induced at OD₆₀₀ ~0.3–0.4 by addition of 0.1 mM IPTG.

In practice, due to experimental limitations, the data must be cleaned. We set any ratios wherein the denominator was too high to measure to zero, assuming that such molecules will not contain substantial desirable fragments. Due to the 2-fold error in the MIC measurements, it was possible to compute ratios greater than one, which were set to one to avoid inclusion of erroneous information. Finally, certain measurements were given as ranges of MIC values and, in each of those cases, we chose the mean of the range.

2.3. Feature selection

The feature selection portion of the Hunting FOX algorithm consists of two steps: (i) permissive LASSO regularization to eliminate non-predictive variables and (ii) hierarchical cleansing to eliminate remaining redundancies. First, we split the data into five disjoint subsets and fit a sparse multinomial classifier with the LASSO penalty⁴⁴ five times, each time reserving one of the five subsets for testing (5-fold cross-validation). We choose to employ multinomial classification due to the natural spread of the data (cf. Fig. S3), and we employ the LASSO penalty as a simple, well-verified way of discarding non-predictive features in the data. Step (i) results in five (overlapping) sets of relevant fragments with associated coefficients from the multinomial fit (see Sec. S1.1) for software details).

From the set of fragments that result in nonzero coefficients in the model, we retain only fragments that have positive coefficients for prediction of class 4 occupancy-that is, fragments whose existence imply a compound will have a high permeating ability-or negative coefficients for prediction of class 0 occupancy-that is, fragments whose existence will decrease the probability of being in the least-permeating class and therefore may be predictive of non-zero permeation. Due to class 0’s higher heterogeneity, no fragments were found with positive coefficients for prediction of class 0 occupancy, so we are not able to present a set of fragments to avoid in addition to a set of fragments to include. We do not retain fragments with coefficients solely pertaining to the middle three classes because the interpretation of such fragments is less clear. We further sparsify the retained fragments by finding a subset that are not hierarchically related (step ii).

Our procedure in step ii is as follows: first, consider only fragments with non-zero coefficients in at least one classifier. Then, find the subset $\left\{\mathcal{H}\right\}$ that are hierarchically related within all classifiers, meaning that for any given fragment $x_i\in \left\{\mathcal{H}\right\}$ appearing in one classifier, either that same fragment or a fragment that it contains or that it is contained by appears in every other classifier. Next, consider fragments in order of increasing coefficient magnitude, where we consider the maximum coefficient magnitude of all the coefficients in all classifiers in which the fragment appears. Each time such a fragment is part of a hierarchical relationship with another fragment, remove it from consideration. Continue doing so until all that remains are fragments that are not hierarchically related to one another. This whole procedure preferentially retains fragments that are ranked as important to at least one classifier and also ensures that all retained fragments contain information present in all classifiers.

2.4. Non-sparse regression model

Having identified a set of likely active fragments, we train a new set of non-sparse classifiers, using only the non-hierarchically-related subset, that may be used to identify molecules from an external library. We employ the same train/test split as before to avoid cross-contamination and fit five multinomial classifiers with a ridge penalty. We run this set of five classifiers on an external testing set; any molecule that is predicted to have an MIC ratio $\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}>0.8$ by all five classifiers is returned as a hit by this run of the Hunting FOX algorithm. We employ five different classifiers to leverage the power of ensemble methods⁴⁵.

2.5. Coarse-grained molecular dynamics

To study the possible molecular contributions of the fragments, we performed biased coarse-grained molecular dynamics simulations employing the MARTINI force field to calculate the two-dimensional potentials of mean force of five molecules crossing a membrane mimic. The molecules we studied were amoxicillin, difloxacin, sarafloxacin, and two molecules from our curated database we refer to as OU-315 and OU-314 (cf. Fig. 4a–b). The parametrization of all molecules studied follows the general MARTINI philosophy, which consists of reproducing the partition coefficient between an organic solvent and water. In our case, the parameters of all drugs were tuned in order to match predicted octanol-water LogP values. (We have attached the itp files as part of the Supporting Material in the folder anc/TOPOLOGIES). In order to retain internal dynamics at the coarse-grained (CG) level, we have employed atomistic generated data using the General Amber force field (GAFF) with charges obtained using the RESP approach⁴⁶. Such simulations were later used to incorporate the necessary bonded terms into the CG geometry as previously performed for similar molecules⁴⁷.

The simulated system consisted of a charge-neutralized outer membrane patch solvated with the polarizable water version of the MARTINI force field. We represented the OM using inhouse developed parameters⁴⁸, based on the GROMOS 53a6 GLYC. The outer leaflet is composed of a homogenous mixture of pure LPS derived from the P. aeruginosa PAO1 BandA strain, which was neutralized using CA++ cations. The inner leaflet is composed purely of DPPE lipids. Although it is possible in general for other lipids to be present in small amounts, this particular composition allows rapid equilibration and has previously been shown to be a good mimic of the P. aeruginosa outer membrane⁴⁹. The system was simulated under constant ionic concentration of 150mM Na-Cl. The final membrane consisted of 12 LPS molecules and 36 DPPE lipids in the inner leaflet. Each compound was placed separately in the water-membrane interface and equilibrated at 310 K under semi-isotropic pressure coupling using a velocity-rescaling thermostat⁵⁰ and a Berendsen barostat⁵¹ respectively for 1 μs before production simulations.

All simulations were carried out with GROMACS 5.1.2⁵², compiled with the PLUMED v 2.5 software⁵³ in order to compute the membrane translocation free energy. Simulations used a 25 fs time-step. Particle Mesh Ewald electrostatics were used with a Coulomb cut-off of 1.1 nm and dielectric constant adjusted to ϵ = 2.5 in order to maintain consistency with the MARTINI polarizable water model⁵⁴. A cut-off of 1.1 nm was used for calculating Lennard Jones interactions, and the Van der Waals potential was shifted to zero at the cut-off. The OM leaflets were coupled to a constant thermal bath maintained at 310 K using a velocity-rescaling thermostat⁵⁰ with a relaxation time of 1.0 ps. Compounds and solvents were coupled separately, due to the needs of the biasing algorithm (see below).

The free energy for membrane translocation was computed using parallel well-tempered metadynamics (PtWtMET)^55,56 (see Sec. S3 for details), requiring a set of four coupled different temperatures. For each drug, four simulations running at 310, 410, 610 and 1010 K were necessary in order to properly converge and overcome the energetic barriers. Production runs were adjusted to one of the coupling temperatures using a NVT ensemble. In order to prevent membrane disruption at higher temperatures, the phosphates of both the Lipid-A region and the DPPE inner leaflet lipids were position-restrained along the δ vector after equilibration.

In our free energy calculations, we chose to bias two specific reaction coordinates: (i) the center of mass (COM) of the drug with respect to the COM of the membrane and (ii) an angle defining the relative orientation of the drug with respect to the membrane (see Fig. S4). Specifically, in the case of amoxicillin, we bias the angle formed between the carbonyl group, the hydroxyl group in the aromatic ring and the COM of the membrane. For both fluoroquinolones, we bias the angle formed between the carbonyl, the distal nitrogen in the piperazine group and the COM of the membrane. Lastly, for OU-315 and OU-314, we bias the angle formed by the nitrogen in the leucoline ring, the central nitrogen in the distal di-amine groups and the COM of the membrane. In all cases, sigma values of 0.1 nm and 0.35 rads were applied to the COM-COM distance and the angle respectively. We find that a bias factor of k = 15 kJ/mol was enough to converge the simulations. We find that a combination of a CG model with the PtWtMET approach allows good convergence at ranges between 12–15 μs of MARTINI time scales.

3.1. Hunting FOX discovers active submolecular fragments

We run the Hunting FOX algorithm twenty-eight times starting from a different random seed each time, which corresponds to a different disjoint split for the 5-fold testing data, and each time identify a set of important molecular fragments (our “chemical vocabulary”), which we employ to train a second set of non-sparse classifiers. We believe this number of repetitions is sufficient to observe the statistics of the algorithm’s behavior. In Fig. 3, we report the number of fragments out of forty-eight that appear between one and twenty-eight times, while in Supporting File anc/top_fragments.xlsx we report all identified fragments and their frequency of appearance. We note a decent level of robustness in the algorithm: although no fragments are reported by every run, one out of forty-eight does appear in twenty-six runs, and twenty-nine fragments appear in more than one run. Although the lack of more fragments appearing repeatedly indicates some dependence on the training data, considering every fragment separately allows for a clear interpretation of the importance of each fragment individually. It might be possible to ameliorate this dependence by employing a representation that considers fragment similarity explicitly (cf. Sec. S4 in the Supporting Information); however such work is outside the scope of the current study.

There are several notable chemical features of the fragments that provide a posteriori empirical support for our procedure. First, twenty-six of forty-eight of all fragments and six of nine fragments appearing more than five times contain part of a benzene ring or one or more whole benzene rings. Such fragments have recently been demonstrated to improve permeability¹², and this result also highlights that the predictive algorithm is able to discriminate 3-dimensional features that are not directly employed in our code. In fact, the presence of rigid benzene rings dramatically improves both rigidity and globularity, important for membrane translocation¹². Second, a large number of the fragments identified contain primary and secondary amine groups. Such groups are expected to strongly interact with both the 2-Keto-3-deoxy-octonate (KDO) and lipid-A regions of the lipopolysaccharide (LPS), functioning as specific anchors for highly anionic membranes (e.g. bacterial OM)⁵⁷. Indeed, there is external evidence that including such groups does improve intracellular accumulation¹² and specificity for Gram-negative bacteria⁵⁸. Third, we note the appearance of a trifluoromethyl fragment (Fig. 3, molecule at a_freq = 9) in nine out of twenty-eight runs. Nowadays, fluorine containing compounds are synthesized in pharmaceutical research on a routine basis and about 10 percent of all marketed drugs contain a fluorine atom. The major rationale is that the presence of fluorine atoms in biologically active molecules can enhance their lipophilicity and thus their uptake and transport. In particular, the trifluoromethyl group (−CF3) confers increased stability and lipophilicity in addition to its high electronegativity. This enrichment with trifluoromethyl-containing fragments is somewhat to be expected, as our initial dataset contained a not-insubstantial number of compounds with this functional group. This provides an additional check of our algorithmic approach, and also points to a potential avenue for extension. It is remarkable that we were able to capture such behaviors without leveraging chemical expertise in the initial steps.

3.2. Molecular Dynamics simulations provide mechanistic rationale of fragment identification

We choose two representative molecules (OU-315, with an experimental MIC ratio of $\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}~0.25-0.5$, and OU-314, with an experimental MIC ratio of $\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}~1.0$) from the training set that summarize the different identified relevant fragments and probe their mechanism of permeation in a relatively high-resolution manner to investigate possible molecular contributions of the fragments selected by the Hunting FOX algorithm in more than five iterations. In Fig. 4a–b we show the candidate molecules as well as highlighted fragments from the set of nine that were identified by the Hunting FOX algorithm. By using a coarse-grained two-dimensional metadynamics calculation in a previously-described membrane model⁵⁹, we compute a potential of mean force (PMF) that suggests mechanistic details by which the fragments enhance the permeation of the selected candidate.

In addition, we choose three previously studied molecules with a range of MIC ratios for comparison: amoxicillin (MIC ratio $\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}~0.015625$), difloxacin, and sarafloxacin (MIC ratios $\frac{{\mu }_{P\Delta 6-\text{Pore}}}{{\mu }_{P\Delta 6}}~0.25$)¹⁴. Of these molecules, difloxacin contains a single one of the top nine fragments, whereas amoxicillin and sarafloxacin contain none. Amoxicillin exemplifies a very polar drug, while both difloxacin and sarafloxacin are fluoroquinolones containing aromatic rings in combination with halogen groups.

A one-dimensional projection of the permeation free energy of the studied compounds is provided in Fig. 4e. Comparison of OU-315 and OU-314 to the other three compounds illustrates two noteworthy features: i) a favorable energy basin in both the core region of the LPS and in the phosphate groups in the 1,2-Dipalmitoyl-sn-glycero-3-phosphoethanolamine (DPPE) lipids of the inner leaflet, and ii) the absence of an energy barrier within the hydrophobic region of the membrane: indeed, the PMF is approximately flat all the way from the outer leaflet to the inner leaflet. The first feature leads to the conclusion that the two compounds containing many favorable fragments are attracted by the anionic regions of the membrane, providing a direct advantage in terms of partitioning within the LPS leaflet. The second feature, the lack of a second energy well within the hydrophobic core of the membrane compared to the three compounds not containing many favorable fragments, provides a potential advantage in terms of permeation by reducing the number of barriers the drug must navigate. We note in addition the lack of these features in three non-permeating control molecules containing few or none of the chemical vocabulary.

The traced approximate lowest free energy paths (red lines) of OU-315 and OU-314 on two-dimensional PMFs including orientation as the second variable (Fig. 4c–d), demonstrate weak zigzag patterns, which suggests that these compounds bypass the aliphatic region via an orientational rearrangement, somewhat resembling the flip-flop mechanism of lipids, and that this mechanism occurs, to within thermal fluctuations, at constant free energy. A plausible explanation for this observed property in the PMF is supported by visual inspection of the simulation of OU-315 (Fig. 4f), which shows water molecules forming a solvation shell within the proximity of the drug, which is expected to overall reduce the penalty of translocating charged groups across the aliphatic region.

3.3. Theoretical and experimental verification of Hunting FOX predictions

We now validate our algorithm theoretically and experimentally, through showing that it is capable of producing well-performing classifiers that can identify novel molecules with desired properties from a library of molecules with unknown properties using only its identified predictive fragments. We train a set of non-sparse classifiers (cf. Sec. 2.4) over 28 iterations of the Hunting FOX algorithm and employ them to make predictions on an external library.

3.4. Conclusions

In this article, we provide proof of concept of the Hunting FOX algorithm, which combines traditional machine learning approaches with a chemical vocabulary-based molecular description inspired by natural language n-grams and FBDD. Despite employing no a priori expert input, it identifies the chemical submolecules that impart compounds with desirable properties and identifies existing drugs with those properties from an external library whose size makes it experimentally intractable, leading to a unique avenue for rational hybrid drug design and drug reuse. Specifically, we have employed our algorithm to identify a set of fragments expected to confer on drugs the ability to permeate the OM of P. aeruginosa, as well as nine compounds expected to be good OM permeators, of which thus far five have been directly experimentally validated. We have also used biased MD simulations to determine the mechanism of permeation of two molecules containing many of the top reported fragments, thus uncovering their molecular-level importance. Hunting FOX opens new portions of chemical space through new combinations of fragments, and also serves as useful tool to generate novel low-molecular weight fragments governing specific target or activity in fragment-based drug discovery. As demonstrated here, this work represents an important step forward for rational hybrid drug design, particularly for antibiotics against Gram-negative bacteria. This simple chemical vocabulary-based methodology generalizes easily to any response variable of interest: for example, one might classify drugs based on their ability to inhibit the growth of cancerous cells compared to that of normal cells and thus identify the reusable chemical vocabulary pertaining to preferential tumor-binding in a rational manner. In future, we plan to employ synthetic biology techniques to design unique hybrid antibiotics from a medicinal chemistry perspective based on the fragments identified herein.