Multiple in silico approaches have hitherto been employed to predict lethal genes and to assess growth metrics on different carbon/nitrogen sources for severe pathogenic organisms including Mycobacterium tuberculosis^{44, 45} and Staphylococcus aureus^{46, 47}. In 2013 and 2014, two studies were published that examined the phenome and gene essentialities of the strain ATCC 17978^{48, 49}. We used these datasets to evaluate the overall performance (functionality and accuracy) of iACB23LX.

Our first validation experiment evaluated the accuracy of our model’s carbon and nitrogen catabolism potentials. More specifically, we detected compounds that could serve as sole carbon and nitrogen sources using the large-scale phenotypic data by Farrugia et al.⁴⁸. Although the authors tested more compounds in total, we could examine only 80 compounds as carbon sources and 48 metabolites as sole nitrogen sources. For the remaining molecules, either no Biochemical, Genetical, and Genomical (BiGG) identifier could be found, or they were not part of the metabolic network. According to the experimental protocol followed by Farrugia et al., we applied the M9 medium and enabled D-xylose as the sole carbon source for the nitrogen testings. As D-xylose was initially not part of the reconstructed network, we conducted extensive literature and database search to include missing reactions. This improved the prediction accuracy, especially for the carbon sources, where an amelioration of 19 % was achieved. In more detail, despite the comprehensive manual curation, the first draft model was reconstructed using the automated tool CarveMe³¹. This resulted in the incorrect inclusion of transport reactions, which were consequently removed to reduce the number of false positive predictions. In both cases, our main objective was to improve the accuracy while keeping the number of orphan and dead-end metabolites low and removing only reactions with no gene evidence (lack of assigned GPR). Similarly, missing reactions were identified and included in the network to eliminate the false negative predictions. For instance, in accordance with the phenotypic data, the strain ATCC 17978 should not be able to grow when utilizing D-trehalose as the sole carbon source. Our model initially predicted a growth phenotype for this carbon source. To overcome this conflict, we deleted the reaction TREP with no organism-specific gene evidence, meaning no GPR was assigned. However, it was not feasible to resolve all inconsistencies since adding transport reactions to resolve false positives or false negatives in the nitrogen testings led to more false predictions in the carbon sources. More specifically, when adenosine, inosine, L-homoserine, and uridine are utilised as sole carbon sources, the model should not predict growth, while sole nitrogen sources should result in a nonzero objective value. In this case, adding transporters would resolve false predictions in the nitrogen tests, while it would have induced more false predictions in the carbon sources tests. All in all, iACB23LX exhibited an overall accuracy of 86.3 % for the carbon and 79.2 % for the nitrogen sources test (Figure 4c); however, after further curation, the accuracy was remarkably improved and reached 87.5 %. By adding their corresponding transport reactions, we resolved discrepancies regarding uridine, inosine, adenosine, and L-homoserine.

Our model was able to catabolise 49 sole carbon and 40 sole nitrogen sources (see Figure 4 (a) and Figure 4 (b), recapitulating totally 69 and 38 experimentally-derived phenotypes, respectively.

We further assessed the ability of our model to predict known gene essentialities. First, 1,164 in silico single gene deletions were conducted on both LB and rich growth media, respectively, to identify all lethal gene deletions. Subsequently, the ratio between the growth rate after and before the respective knockouts (FC_gr) was calculated, and the genes were accordingly classified (see Materials and Methods). For the optimization, two mathematics-based approaches from the Constraints-Based Reconstruction and Analysis for Python (COBRApy)⁵⁰ package were deployed: the FBA⁵¹ and the Minimization of Metabolic Adjustment (MOMA)⁵². Between the two methods, a similar distribution of the FC_gr values was observed (Figure 5 (a) and Figure 5 (b)). Using FBA, 97, 75, and 991 genes were predicted to be essential, partially essential, and inessential on the LB medium, respectively, whereas optimization with MOMA resulted in 110, 85, and 968 genes (Figure 5 (c) and Supplementary Files S2 and S3). These genes were primarily associated with the biosynthesis of cofactors and vitamins, the amino acid/nucleotide metabolism, the energy metabolism, and the metabolism of terpenoids and polyketides. Additionally, we examined in more detail how nutrition availability impacts the gene essentiality by conducting single-gene knockouts in the rich medium. Both optimization methods resulted in more essential genes when the model was required to alter its metabolic behavior due to the absence of nutrients, i.e., with the LB growth medium, compared to the rich medium (Figure 5 (c) and Supplementary Files S2 and S3). In general, FBA detected more genes to be dispensable for growth in both nutritional environments. On the other hand, MOMA classified more genes as essential or partially essential (Figure 5 (c) and Supplementary Files S2 and S3), while genes from FBA build a subset of the essential genes derived by MOMA. Furthermore, we validated the prediction accuracy of iACB23LX using already existing gene essentiality data. At the time of writing, the transposon mutant library by Wang et al. is the only ATCC 17978-specific experimental dataset⁴⁹. With this dataset together and the LB medium, our model demonstrated an accuracy of 87 % with both optimization methods (Figure 5 (d)). We further analysed the predicted false negative genes and probed their proteomes to investigate the existence of human orthologs (see Materials and Methods and Supplementary Table S4). With this, we aimed to eliminate crosslinkings to human-similar proteins since metabolic pathways or enzymes that are missing from the human host have been an important resource of druggable targets against infectious diseases⁵³. From the 37 genes that our model predicted to be essential (with FBA and MOMA) contradicting the experimental results, 17 were found to be non-homologous (see Supplementary Table S4). Some examples are the genes encoding the enolpyruvylshikimate phosphate (EPSP) synthase (A1S_2276), chorismate synthase (A1S_1694), riboflavin synthase (A1S_0223), phosphogluconate dehydratase (A1S_0483), and 2-keto-3-deoxy-6-phosphogluconate (KDPG) aldolase (A1S_0484). The EPSP synthase converts the shikimate3-phosphate together with phosphoenolpyruvate to 5-O-(1carboxyvinyl)-3-phosphoshikimic acid. Subsequently, the chorismate synthase catalyses the conversion of the 5-O-(1carboxyvinyl)-3-phosphoshikimic acid to chorismate, the seventh and last step within the shikimate pathway⁵⁴. Chorismate is the common precursor in the production of the aromatic compounds tryptophan, phenylalanine, and tyrosine, as well as folate and menaquinones during the bacterial life cycle. The shikimate pathway is of particular interest due to its absence from the human host metabolome and its vital role in bacterial metabolism and virulence. Moreover, the enzyme riboflavin synthase catalyses the final step of riboflavin (vitamin B2) biosynthesis with no participating cofactors. Riboflavin can be produced by most microorganisms compared to humans, who have to externally uptake them via food supplements. Also, it plays an important role in the growth of different microbes, especially due to its photosynthesizing property that marks it as a non-invasive and safe therapeutic strategy against bacterial infections⁵⁵. Lastly, the phosphogluconate dehydratase catalyses the dehydration of 6-Phospho-D-gluconate to KDPG, the precursor of pyruvate and 3-Phospho-D-glycerate⁵⁶. This enzyme is part of the Entner–Doudoroff pathway that catabolises glucose to pyruvate, similarly to glycolysis, but using a different set of enzymes⁵⁷.

We further assessed the druggability of our essential nonhomologous proteins and investigated the existence of inhibitors or compounds known to interact with the enzymes. For this, we used the online DrugBank database that contains detailed information on various drugs and drug targets⁵⁸. In all cases, the listed drugs are of unknown pharmacological action, and there is still no evidence indicating the enzymes’ association with the molecule’s mechanism of action. For instance, the flavin mononucleotide and the cobalt hexamine ion were listed as known inhibitors of yet unknown function against the chorismate synthase, while glyphosate, shikimate3-phosphate, and formic acid have been experimentally found to act with EPSP synthase. Six non-homologous genes were marked as hypothetical or putative in the KEGG³⁰ database and/or lacked enzyme-associated information. We searched for drug leads by aligning the query sequences against the DrugBank’s database to find homologous proteins. Two out of six were found to have a protein hit. More specifically, the protein encoded by A1S_0589 was found to have high sequence identity with the phosphocarrier protein HPr of Enterococcus faecalis (Bit-score: 48.5), while the translation product of A1S_0706 resembles the sugar phosphatase YbiV of Escherichia coli (Bit-score: 225.3). According to DrugBank, dexfosfoserine and aspartate beryllium trifluoride have been experimentally determined to bind to these enzymes; however, their pharmacological action is still unknown. The Supplementary Table S4 lists all non-homologous essential genes reported for iACB23LX.

Overall, iACB23LX exhibits high agreement to all validation tests and can, therefore, be used to systematically derive associations between genotypes and phenotypes.