Morphological profiling of tubercule bacilli identifies drug pathways of action

Morphological profiling is a method to classify target pathways of antibacterials based on how bacteria respond to treatment through changes to cellular shape and spatial organization. Here, we utilized the cell-to-cell variation in morphological features of Mycobacterium tuberculosis bacilli to develop a rapid profiling platform called Morphological Evaluation and Understanding of Stress (MorphEUS). MorphEUS classified 94% of tested drugs correctly into broad categories according to modes of action previously identified in the literature. In the other 6%, MorphEUS pointed to key off-target or secondary bactericidal activities. We observed cell-wall damaging activity induced by bedaquiline and moxifloxacin through secondary effects downstream from their main target pathways. We implemented MorphEUS to correctly classify three compounds in a blinded study and identified an off-target effect for one compound that was not readily apparent in previous studies. We anticipate that the ability of MorphEUS to rapidly identify pathways of drug action and the proximal cause of bactericidal activity in tubercule bacilli will make it applicable to other pathogens and cell types where morphological responses are subtle and heterogeneous. Significance Statement Tuberculosis is a leading cause of death in the world and requires treatment with an arduous multidrug regimen. Many new tuberculosis drugs are in development, and the drug development pipeline would benefit from more rapid methods to learn drug mechanism of action and off-target effects. Here, we describe a high throughput imaging method for rapidly classifying drugs into categories based on the primary and secondary cellular damage called Morphological Evaluation and Understanding of drug-Stress (MorphEUS). We anticipate that MorphEUS will assist in rapidly pinpointing causes of cellular death in response to drug treatment in tuberculosis and other bacterial pathogens.


71
Mycobacterium tuberculosis (Mtb), the causative agent of tuberculosis (TB), causes more deaths 72 annually than any other infectious agent (1). Tuberculosis treatment is lengthy, lasting between four months 73 to over a year (1). The difficult regimen, rate of relapse, and incidence of drug resistant Mtb has motivated 74 a significant effort to develop new antibacterial compounds that are effective in sterilizing Mtb infection (2).

75
Many new drug classes and derivative compounds have been developed (2), but rapidly identifying the 76 primary and secondary pathways of action of each compound is often a protracted process due to the 77 difficulty in generating resistant mutants, and dissecting the broad-reaching metabolic effects of drug 78 treatment leading to death (3). Furthermore, drug action on bacterial cells can elicit dynamic responses in 79 multiple pathways both on-and off-target, some of which are specific to bacterial growth environment and 80 treatment dose (4-6). An approach that combines computational algorithms with traditional methods for 81 identifying pathway of action would enable faster drug development.

82
In other bacterial systems such as Escherichia coli, Bacillus subtilis, and Acinetobacter baumannii, 83 profiling of cytological changes in response to treatment has yielded a rapid and resource-sparing 84 procedure to determine drug mechanism (7-9). This method, known as bacterial cytological profiling (BCP), 85 groups drugs with similar mechanisms of action by clustering profiles of drug-treated bacteria using 86 multivariate analysis methods such as principal component analysis (PCA) (7-9). BCP is efficient and rapid 87 because these cytological features can be automatically derived from high-throughput images of stained, 88 fixed samples.

89
We hypothesized that BCP could be similarly utilized to map pathways of drug action in Mtb. We 90 found that Mtb morphological response to treatment were subtle and incorporation of cellular variation 91 metrics improved the ability to profile drug action in Mtb. Here, we present an imaging and data analysis 92 pipeline for Mtb that classifies morphological profiles called MorphEUS (Morphological Evaluation and 93 Understanding of drug-Stress). We demonstrate that MorphEUS clusters antibacterials by their pathways 94 of action. Because MorphEUS is based on the physical unraveling of cells due to stress, it can be used to 95 identify key proximal cellular stressors that arise from off-target and secondary effects. We used MorphEUS 96 to identify secondary effects for two TB drugs (moxifloxacin and bedaquiline) and a non-commercial 97 4 compound. We propose that MorphEUS will be useful in classifying drug action in other pathogens where 98 morphological responses are subtle and heterogeneous.

100
Antibacterial treatment induces drug-specific morphological response in mycobacteria

101
To assess whether mycobacteria exhibit distinct morphological responses to drugs targeting 102 different pathways, we used a chromosome-decorating reporter strain (GFP translational fusion to RpoB) 103 of Mycobacterium smegmatis to visualize cell and nucleoid shape characteristics (10) before and during 104 drug treatment by time-lapse imaging. We observed rapid nucleoid condensation in rifampicin-treated M.

105
smegmatis that was not apparent in ethambutol-treated cells (Fig.1A). In contrast, moxifloxacin treatment 106 induced filamentation with nucleoid decompression. These antibacterials have different mechanisms of 107 action (inhibitors of transcription, cell-wall synthesis, and DNA synthesis, respectively), providing evidence 108 that the morphological changes induced by drug treatment in mycobacteria are dependent on drug target.

109
These data led us to hypothesize that exploitation of morphological features in mycobacteria would allow 110 clustering of antibacterials in a manner similar to previously described BCP studies.

111
We next asked whether drug treatment of Mtb, like M. smegmatis, elicited well-defined cytological 112 fingerprints. Guided by cytological profiling methods in other bacterial species (7-9), we treated Mtb grown 113 in standard rich growth medium with a high drug dose (3x the 90% inhibitory concentration; IC90) for 17 114 hours (~1 doubling time). We generated a dataset of morphological features from Mtb treated with 34 115 antibacterials that encompass a wide range of drug classes and classified the target pathway of each drug 116 according to published findings (Table S1). We imaged fixed, membrane-(FM4-64FX) and nucleoid-(SYTO  Table S1). Furthermore, the morphological responses to drug treatment in Mtb were more subtle than M. 120 smegmatis (Fig. 1A). Because Mtb are smaller than M. smegmatis, we hypothesized that there may still be 121 quantitatively significant differences among morphological features in drug-treated bacilli that are more 122 sensitive to noise and variation. Using image segmentation and analysis, we quantified 25 morphological 123 features per treatment group (Table S2). We observed significant differences among treatment groups in 124 features such as cell shape, nucleoid shape, and staining intensity (Fig. 1B). However, the resulting 125 5 morphological profiles from drug treatment did not cluster based on broad drug target categories using PCA 126 ( Fig. S1A top).

127
One explanation for the poor performance of BCP in Mtb may be the significant cell-to-cell variation 128 in morphological features ( Fig. 1 and Fig. S2). This inherent heterogeneity is consistent with the variable 129 nature of Mtb, which on the single-cell level exhibits heterogeneity through asymmetric growth and division, 130 differential drug susceptibility, and metabolic state (11)(12)(13)(14)(15). Cell-to-cell variation is most apparent in the 131 ability of Mtb bacilli to take up stains. For example, only ~10% of untreated bacilli are stain positive, whereas 132 the proportion increases to ~30% when treated with cell-wall acting antibacterials (Fig. 1B). We speculated 133 that variation itself was an important feature of drug response that should be captured in the profiling of 134 drug mechanism of action. To capture cell-to-cell variation, we developed a new analysis pipeline for Mtb that incorporates 138 variation as an important class of features to discriminate drug target pathways (Table S1- (16)). This analysis formulation also addresses the subtlety of cytological changes by taking into 140 account the full dimensionality of the data to produce discrete classifications. The exploitation of feature 141 variation provided another dimension to distinguish drug categories. For example, when treated with 142 isoniazid, Mtb nucleoid stain intensity was less variable than when treated with bedaquiline or meropenem 143 ( Fig. 1B middle). We accounted for a fragile feature selection process (in which several sets of features 144 may achieve similar model accuracy) by performing a series of classification trials (Fig. 2). The resulting 145 analysis was visualized using a network web or matrix describing the frequency of drug-drug links (Fig. 3).

146
Morphological responses were measured in cells treated with high-and low-doses of drugs (Fig. S4) and 147 combined to generate a joint-dose profile to incorporate richer profiling behaviors. We refer to this analysis 148 pipeline as MorphEUS (Morphological Evaluation and Understanding of drug-Stress).

149
Using MorphEUS, drugs in the same broad categories are generally grouped together (94% 150 accurate categorization with 76% accurate cross validation; Fig. 3). Furthermore, within the broad 151 categories, drugs sharing target pathways were found to have stronger connections to one another. For 152 example, within the cell-wall acting category, strong connections were observed between ethionamide and 6 isoniazid (inhibitors of the enzyme InhA), delamanid and pretomanid (nitroimidazole drug class, inhibits the 154 synthesis of mycolic acids), and meropenem, cefotaxime, and vancomycin (all peptidoglycan inhibitors) 155 ( Fig. 3 and Table S1). We also observed strong connectivity between inhibitors of cellular respiration with 156 the ionophores CCCP, monensin, and nigericin forming stronger connections with each other compared to 157 clofazimine and thioridazine (both shown to target NDH-2 of the electron transport chain) (17, 18).

158
MorphEUS also predicted strong connections among protein synthesis inhibitors that target the 50S 159 ribosomal subunit (clarithromycin, chloramphenicol, and linezolid). Finally, the fluoroquinolones levofloxacin 160 and ofloxacin grouped together as did rifampicin and rifapentine -inhibitors of transcription.  S6). We focused on induction of the iniBAC operon, which is rapidly upregulated due 171 inhibition of cell-wall synthesis and is used to screen for cell-wall acting compounds (12, 13). We observed 172 that cycloserine's induction of iniBAC genes was mild compared to other cell-wall acting antibiotics (~1.5 173 fold and ~3.5 fold respectively), consistent with the weak association of cycloserine to other cell-wall 174 antibacterials ( Figure 3 Right and Fig S6 Top). We next hypothesized that cycloserine has off-target effects 175 that are DNA damaging and that these effects drive the connection to moxifloxacin. However, we did not 176 find reports of mutants that confer resistance to both moxifloxacin and cycloserine (14-18). An alternative 177 hypothesis is that the connection cycloserine to moxifloxacin is driven by off-target effects by moxifloxacin, 178 such as cell-wall damage. Moxifloxacin's morphological profiles are highly dose-dependent, resembling the 179 other fluoroquinolones at low-dose but not at high-dose or joint-dose (Fig. S7, Fig 3). This shift away from 180 other fluoroquinolones by moxifloxacin suggests that there is an off-target or secondary effect at high dose.

7
Transcriptional analysis based on published studies of moxifloxacin-treatment Mtb demonstrated a mild but 182 significant increase in iniB expression and a 2.5-fold induction in ddlA, a molecular target for cycloserine 183 (Fig. S6 bottom) (16,17,19). Taken together these data suggest that the similarity between cycloserine 184 and moxifloxacin morphological profiles arises from an off-target cell-wall damaging effect of moxifloxacin 185 and a mild inhibition of cell-wall synthesis by cycloserine.

186
The second unexpected profile was for bedaquiline, an ATP-synthesis inhibitor, which mapped to

211
Blinded to compound identity, we next used MorphEUS to identify pathways of action for three non-212 commercial antituberculars with known mechanisms of action. We mapped unknown compounds 1 and 2 213 as cell-wall acting; compounds 1 and 2 were nearest neighbors to ethionamide and ethambutol, respectively 214 ( Fig. 4). We unblinded the compound identities to compare their known mechanisms of action to those 215 predicted by MorphEUS. These compounds (DG167 and its derivative JSF-3285) were validated through 216 extensive biophysical, X-ray crystallographic, biochemical binding and spontaneous drug-resistant mutant 217 studies to be inhibitors of cell-wall mycolate biosynthesis through specific engagement of the a-ketoacyl 218 synthase KasA ( (27), and Cell Chem. Biol, provisionally accepted). Taken together, our analysis of DG167 219 and JSF-3285 using MorphEUS has independently validated the target pathway of these two compounds 220 and shown these analogs act through the same pathway of action.

221
The activity of the third unknown compound was harder to interpret. MorphEUS predicted unknown 222 compound 3 to act on both the cell-wall and DNA and had ofloxacin as its nearest neighbor (via joint-dose 223 profiles; Fig. 4). In contrast, MorphEUS analysis at low treatment dose mapped unknown compound 3 to 224 cell-wall acting antibacterials with pretomanid as its nearest neighbor. Unknown compound 3 therefore 225 displays dose dependent effects that become apparent in the joint profile suggesting that downstream off 226 target effects are amplified upon increasing the treatment dose. We unblinded the compound to learn if our 227 conclusions were corroborated with previous mechanistic studies performed. Unknown compound 3 was 228 the recently published triazine JSF-2019 that resembles pretomanid in both its F420-dependent production 229 of NO• and its ability to inhibit mycolic acid synthesis, albeit at a different step in the pathway (28). The

258
In applying MorphEUS to a set of 34 known antibiotics and three blinded, non-commercial 259 antibacterials, we found that MorphEUS grouped antibacterials by their killing activity, which may be the 260 primary target pathway or off-target effects. Most of the antibacterials profiled with their known pathway of 261 action, but we identified three drugs (moxifloxacin, bedaquiline, and JSF-2019) where killing activity 262 mapped instead to off-target effects. Bedaquiline and moxifloxacin are known as inhibitors of respiration 263 and DNA synthesis, respectively, yet both clustered with cell-wall-inhibiting drugs. For bedaquiline, we 264 found that the cell-wall damage effect was specific to metabolic state and was a secondary effect for cells 265 10 carrying out glycolysis and was not observed for cells utilizing fatty acids as a carbon source. We

287
We have shown that the MorphEUS pipeline identifies drug pathways of action, also revealing off-288 target and downstream drug effects that are proximal to antibacterial action. We anticipate that the    were nearest neighbors with each other. We found no likeness (e.g. 100% confusion) between the same 372 treatment groups (e.g. that the controls were not identifiable into similar treatment groups), suggesting that 373 solvents alone did not induce strong morphological effects.

15
Montage images were generated using a custom macro that captures 25 individual fields of view per image.

406
Two technical replicate images were taken from each sample for a total of 50 images per biological 407 replicate. Three biological replicates were generated for each drug treatment. Images were recorded with 408 a DV Elite CMOS camera for all three channels.

426
Transformations were recovered in 1 ml of 7H9-rich-ADC medium, incubated for three hours then plated 427 on 7H10-ADC plates containing hygromycin B at 50 µg/ml. The presence of the C-terminal GFP 428 translational fusion to RpoB was validated by fluorescence microscopy using the FITC (Ex. 475nm Em. 429 525nm) channel as described above. The rpoB targeting oligo sequence was:

503
This transformation is then applied to the entire dataset, including treated and untreated samples (Fig. S3).

504
To perform TVN, we dedicated 25% of our samples in every experiment and imaging session to be 505 untreated controls. Each classification trial begins with stochastically selecting 80 untreated controls; thus, 506 feature sets were restricted to a maximum of 79 features since a PCA transform with n samples can only 507 have n-1 features.

509
Principal Component Analysis. PCA was performed on the normalized data (feature selected or not, as 510 indicated) using the built-in MATLAB function. In the BCP pipeline, PCA reduces the dimensionality of the 511 data, allowing variance across all of the features to be visualized in 3 or fewer dimensions. After accounting 512 for heterogeneity with batch normalization and including features of variation into the profiles, some drug 513 clustered observed, especially among cell-wall acting antibacterials (Fig. S1A lower). By PCA, bedaquiline 514 clustered with cell-wall acting drugs in standard, rich growth medium (Fig. S1A lower and Fig. S5B left) but 515 not in fatty acid-rich growth medium (Fig. S5B right).

516
19 517 K-nearest neighbors. K-nearest neighbors (KNN) analysis was implemented using the cosine distance 518 metric and the knnsearch MATLAB function. K was set to 1, thus only the first nearest neighbor was 519 identified. For our setup, we took the median PCA score from the three replicates for each drug as inputs 520 for the KNN analysis. The KNN algorithm finds the k-nearest neighboring points where the cosine distance 521 between PCA scores is shortest. MATLAB defines the cosine distance as one minus the cosine of the 522 included angle between points. We observed that feature selection was dependent on which untreated

531
Iterative feature selection. To reduce overfitting and noise in our 94 variable feature set, we utilized the 532 Minimum Redundancy Maximum Relevance (mRMR) feature selection algorithm (32). Here we customized 533 previously published MATLAB code to perform mRMR feature selection using the mutual information 534 difference scheme (32). Because the algorithm rank orders variables and does not automate the selection 535 of the optimal number of features, we implemented an iterative feature selection method that rewards runs 536 that result in more drug-drug connections with same target pathways (Table S1). Since we begin with 94 537 features but are limited to 79 variables by our TVN analysis, mRMR was used to rank order the top 79 538 features. Starting with the 79 rank-ordered features, we removed each feature individually and performed 539 TVN, PCA and KNN analysis on the remaining feature set. Success of the feature set was quantified by 540 accuracy of the KNN in linking drugs belonging to the same broad category assigned by literature review 541 (see Table S1). The feature set that resulted in greatest model accuracy was selected, and the variable 542 removal process was repeated until maximal prediction performance was reached. On average, these 543 iteratively determined feature sets contained 38 variables. 544 20 545 Consensus KNN (cKNN): The cKNN results compiled from all 70 classification trials were visualized using 546 a network map and heatmap, where edge color and grid square color, respectively, corresponds to how 547 frequently two drug profiles were identified as nearest neighbors. Because our goal was to identify similar 548 treatment profiles, connections between profiles in the cKNN were made undirected and the drug-drug 549 categorization matrices (such as Fig. 3B left) are symmetric. These visuals allow for classification of the 550 target biological pathway of each drug based on the robustness of phenotypic similarities between drug 551 profiles can easily be evaluated. All maps and plots were generated in MATLAB (2019a).

553
Comparison to random. As a comparison to model accuracy by random change, we tested how accurate

554
MorphEUS was when the drug categories were randomly assigned. To do so, the labels for the drugs in 555 the final cKNN were randomly swapped, resulting in 22% accuracy compared to 94% for the joint dose  Table S1.

566
Low dose, high dose, and joint dose profiles. Mtb cytological features are dependent on drug target but 567 also treatment dose and duration (Fig. S6). This raised the possibility that morphological profiles from a low 568 dose of treatment or a joint profile of low (0.25xIC90) and high (3xIC90) dose treatments would improve the 569 accuracy of drug classification using the full drug set and subsequent cross validation. To investigate 570 whether a joint dose profile best describes the variation in the morphological response in Mtb, the full 94 571 feature datasets from both drug doses were concatenated, resulting in 188 total features. We also applied 572 21 MorphEUS to low dose and high dose treatments as separate profiles. We observed high accuracy using 573 each of the dose treatments (low, high, and joint as 97, 91, and 94% respectively), but the joint dose profiles 574 were better cross validated (76%) compared to high (68%) and low (62%) dose MorphEUS. We therefore 575 use joint dose profiling as the default for MorphEUS.

577
Classification of unknown compounds. We apply new compounds to MorphEUS in the same manner 578 as cross validation, only the MorphEUS pipeline is done on the full 34 drug set and the unknown is added 579 to the set for the final KNN during each classification trial. 580 581 Statistical Analysis. We performed the Kruskal-Wallis test to identify drug treatments that induce 582 significantly different morphological features compared to untreated cells in rich medium ( Fig. 1B and S5A).

583
While the MorphEUS pipeline utilizes population-based features, the Kruskal-Wallis test was applied to the 584 features of individual cells (n=1625-3983 for rich, n=1029-6733 for conditions). The Kruskal-Wallis test was 585 applied to each drug and/or environmental condition individually, per feature. In each case the null 586 hypothesis was that the median feature value for Mtb exposed to a specific drug and/or environmental 587 stress was drawn from the same distribution as the median feature value for the untreated controls.