Dual transcriptome based reconstruction of Salmonella-human integrated metabolic network to screen potential drug targets.

Salmonella enterica serovar Typhimurium (S. Typhimurium) is a highly adaptive pathogenic bacteria with a serious public health concern due to its increasing resistance to antibiotics. Therefore, identification of novel drug targets for S. Typhimurium is crucial. Here, we first created a pathogen-host integrated genome-scale metabolic network by combining the metabolic models of human and S. Typhimurium, which we further tailored to the pathogenic state by the integration of dual transcriptome data. The integrated metabolic model enabled simultaneous investigation of metabolic alterations in human cells and S. Typhimurium during infection. Then, we used the tailored pathogen-host integrated genome-scale metabolic network to predict essential genes in the pathogen, which are candidate novel drug targets to inhibit infection. Drug target prioritization procedure was applied to these targets, and pabB was chosen as a putative drug target. It has an essential role in 4-aminobenzoic acid (PABA) synthesis, which is an essential biomolecule for many pathogens. A structure based virtual screening was applied through docking simulations to predict candidate compounds that eliminate S. Typhimurium infection by inhibiting pabB. To our knowledge, this is the first comprehensive study for predicting drug targets and drug like molecules by using pathogen-host integrated genome-scale models, dual RNA-seq data and structure-based virtual screening protocols. This framework will be useful in proposing novel drug targets and drugs for antibiotic-resistant pathogens.

Introduction

S. Typhimurium, is a gram-negative invasive and facultative pathogen that can infect various animal species [1]. Upon infection, it mostly causes food poisoning and leads to gastroenteritis in humans [2]. Salmonella infection causes 130,000 deaths every year, and it mostly affects people in low income countries [3]. S. Typhimurium is an intracellular pathogen, residing inside a membrane-bound compartment within host cells during infection [4]. This compartment is called Salmonella-containing vacuole (SCV), and it enables proliferation of S. Typhimurium inside the host cell by escaping from the host defense mechanism. Although SCV is considered a nutrition poor environment, this does not pose a problem for S. Typhimurium due to its highly adaptive lifestyle [1]. Besides being a highly adaptive pathogen, the increasing rate of antibiotic resistance by S. Typhimurium has the potential to become a serious concern for public health. Therefore, it is important to identify novel drug targets to eliminate S. Typhimurium infections [5]. Understanding the metabolic activities of a pathogen is an important part of drug development process. Transcriptome data is a key data type that can elucidate metabolic activities of an organism [6]. It reflects enzymatic activities of the organism via mRNA levels that carry genetic information of those enzymes. One of the novel approaches in transcriptomics is dual RNA-sequencing (dual RNA-seq), which can be used to elucidate pathogen-host relationship since it measures mRNA levels of pathogen and host simultaneously during infection [7, 8].

Genome-scale metabolic network (GMN) models have shown utility in the analysis of metabolic activities of pathogen and host during infections [7, 9] Analysis of GMN models with constraint-based techniques to predict novel drug targets has the advantages of (i) being cost effective, (ii) being time efficient, (iii) providing a wide range of analyses of metabolic pathways at the same time. There are several studies about the prediction of novel drug targets to eliminate pathogen induced infections by analyzing GMN models [10]. Most widely used constraint-based computational approach for the analysis of GMNs is Flux Balance Analysis (FBA). FBA is a mathematical optimization technique that uses linear optimization to predict distribution of metabolic fluxes at steady state conditions. It uses an objective function besides constraints to select an optimum point from the flux solution space [11]. In silico gene deletion analysis is another widely used constraint-based analysis technique, which is used to determine potential drug targets [12-18]. Several GMN models were reconstructed so far for different Salmonella strains [19-21]. But these models are generic models, and they do not represent the metabolism of Salmonella inside a host cell. There are several techniques to create condition specific GMN models by mapping transcriptome data on to the generic GMNs. One of the commonly used transcriptome data mapping methods to generate condition specific GMNs is Gene Inactivity Moderated by Metabolism and Expression (GIMME) algorithm [22], which predicts active and inactive reactions in a GMN based on mRNA levels belonging to a particular condition.

Infection leads to a set of intricate interactions between pathogen and host cells, and these interactions should be taken into account in the process of identification of novel drug targets [9]. Pathogen-host integrated GMNs have potential to shed light on pathogen-host interactions (PHI) when integrated with dual RNA-seq data [7]. Pathogen-host metabolic modeling is a multi-cellular interaction modelling approach, where GMNs of both pathogen and host organisms are integrated in the simulations of metabolic phenotypes. Even if not commonly used yet, some pathogen-host GMNs are available in the literature [23, 24]. Here, we aim to provide a better insight into S. Typhimurium and host interactions by taking advantage of pathogen-host integrated GMNs and dual RNA-seq data in order to determine novel drug targets that can eliminate S. Typhimurium induced infections. We further report potential drugs for the identified drug target candidates by using a drug repositioning based approach. Through a prioritization criteria, pabB was selected as a high-ranked putative target, and a structure-based screening of novel drugs for pabB was performed using molecular docking simulations. To our knowledge, this is the first study that reconstructs a condition-specific pathogen-host GMN model by mapping dual RNA-seq data to predict novel drug targets.

Results

Pathogen-host integrated genome scale metabolic network analysis

A pathogen-host integrated GMN model was reconstructed in this study for the first time in literature for S. Typhimurium, with a total of 3586 genes from both organisms and 11,029 reactions. This model was used to generate condition-specific pathogen-host GMN models by mapping the dual RNA-seq data from 0th, 8th and 16th hours of infection [25]. GIMME was used to integrate the model and the transcriptome data, and the number of reactions for the corresponding condition-specific pathogen-host GMNs are given in Table 1.

Table 1

The reaction and metabolite numbers of condition specific pathogen-host GMN models.

	Condition-specific GMN at the beginning of infection	Post-infection condition specific GMN at 8th hour	Post-infection condition specific GMN at 16th hour
Number of Reactions	8773	8933	9089
Number of Metabolites	6941	6982	6979
Number of HeLa reactions	6595	6681	6536
Number of S. Typhimurium Reactions	2178	2252	2553

Pathogen-host GMN was created by combining generic human and S. Typhimurium GMN models, and the number of reactions decreased in the condition-specific pathogen-host GMN models since they represent infected HeLa cell by S. Typhimurium in specific conditions. Nearly 2500 reactions were discarded from the generic Pathogen-Host GMN model since they were not active at the beginning of infection based on mRNA levels. On the other hand, the reaction profile of condition-specific pathogen-host GMN is different between the beginning and the late stages of infection. As infection progresses, the number of HeLa reactions remains almost the same while S. Typhimurium reactions dramatically increase (Table 1). The reaction profiles at the beginning of infection (0th hour) and at late stage of infection (16th hour) were compared. Even if there is not a considerable change in the number of HeLa reactions, the reaction profiles are different between the two conditions (Fig 1). During the infection, S. Typhimurium must adapt to the nutrition environment and physical conditions to survive inside the host [26]. Therefore, the increase in the number of S. Typhimurium reactions during the infection can be attributed to the activation of genes that might be necessary for the adaptation and proliferation of the pathogen. The reaction profiles of condition-specific pathogen-host GMNs were compared to each other in order to identify alterations in metabolic pathways based on the progress of infection. For host and pathogen separately, reactions only active in the condition-specific GMN at the beginning of infection (GMN0th) and only active in the post-infection condition-specific GMN at the 16th hour (GMN16th) were identified and grouped by their pathways. The identified metabolic pathways and corresponding number of infection-time specific reactions are given in Fig 1.

Fig 1

The differences in reaction profiles of condition-specific pathogen-host GMNs grouped by pathways.

Pathways with at least 2 differential reactions are given. Red and green bars represent the number of HeLa cell reactions that are only active in 0th hour and 16th hour respectively. Purple and blue bars represent the number of S. Typhimurium reactions that are only active in 0th and 16th hour respectively. The pathway names are listed on the left side of figure.

Fig 1. shows that there is severe modulation in the lipid related pathways of host in infection. Lipid related pathways of the host are known to be subjected to modulation during the invasion of bacteria [27]. On the other hand, multiple metabolic pathways of S. Typhimurium are altered during infection based on Fig 1. Most dramatic change is in glycerophospholipid metabolism, which is one of the most important pathways for dual-membrane envelope of gram-negative bacteria [28]. There is also a dramatic change in alternate carbon metabolism, implying the utilization of different carbon sources other than glucose in the late stage of infection. Even if glucose is major carbon source for S. Typhimurium, it utilizes different carbon sources during infection [29]. Cofactor and prosthetic group biosynthesis, which is tightly related to enzymatic activities, is also altered as expected since enzymatic activities become varied as infection progress.

The first step in the validation of GMN models is comparing predicted flux rates with experimental data from literature. FBA was performed to predict flux rates of GMN models by maximizing TB (see Material and Methods). Ethanol and succinate production rates were set to zero based on literature information [30], and the other constraints were set as detailed in the section Materials and Methods. Predicted secretion rates of major by-products for S. typhimurium are given in Table 2 together with the literature-reported secretion rates at infection. The relative rates of by-products predicted by the pathogen-host integrated condition-specific model at 16th hour are in perfect agreement with the experimental data from HeLa-infecting S. Typhimurium cells [30]. Repeating simulations by using the non-reduced model leads to rates about 20% higher than the reduced model predictions. The use of only reduced Salmonella model without the host network, on the other hand, led to 25% higher acetate secretion rates. Predicted acetate, formate, and lactate secretion rates at 8th hour of infection, on the other hand, are very close to the rates predicted for the beginning of infection (Fig 1). Therefore, for the rest of the study, we used the model reconstructed for the 16th hour of infection as the model representative of the infectious state of the organism.

Table 2

Predicted flux rates obtained by FBA analysis of condition-specific pathogen-host integrated GMNs are compared with the experimental results.

	Infection (0th hr) Flux Values (mmol/gDW/h)	Infection (8th hour) Flux Values (mmol/gDW/h)	Infection (16th hour) Flux Values (mmol/gDW/h)	Experimental (nM/cell/h) [30]
D-Lactate Secretion	0	0	11.29	10 ± 3
Acetate Secretion	1.94	2.17	6.80	4 ± 2
Formate Secretion	10.86	9.14	3.55	2 ± 1
Succinate Secretion	0	0	0	0
Ethanol Secretion	0	0	0	0

Identification of potential drug targets

Prediction of essential genes (EG) for the survival of pathogen inside host organism is the primary step for most of the drug discovery processes. Enzymes produced from EGs are potential drug targets that can be targeted with chemical molecules to eliminate the pathogen. EGs for the infection were predicted by using GMN16th, which represents the infectious state. Here, 140 EGs were predicted for the infection. Predicted 140 EGs were compared with literature by using Database of Essential Genes (DEG) [31], which reports data from three experimental gene deletion studies from rich medium experiments [32-34] Data were available for 137 of 140 predicted EGs, 93 of which were reported as essential genes in at least one study (68%). (S1 Table). The genes falsely predicted as essential can be attributed to the fact that the experiments were performed in rich medium conditions with no host cells involved while the simulations were performed by the pathogen-host integrated genome-scale metabolic network.

Drug targets should not show high amino acid sequence similarity with human proteins in order to prevent side effects. Therefore, homology analysis was performed to identify drug targets that are similar to human proteins. The similarity was determined based on the predefined cutoff value detailed in Materials and Methods. Out of 140 potential drug targets, 52 proteins were discarded since they have high similarity with human proteins (S2 Table). Pathway enrichment analysis was performed with non-homologous 89 proteins in order to characterize potential drug targets (Fig 2). The most enriched pathway is the lipopolysaccharide biosynthesis pathway, which is critical for the survival of S. Typhimurium since it maintains the functionality of the outer membrane of the pathogen [35]. Another enriched pathway is the biosynthesis of amino acids, and it is a reasonable metabolic pathway that can be targeted since it is indispensable for pathogens [36]. Like lipopolysaccharides, peptidoglycans are also indispensable molecules for the functionality of bacterial cell wall [37], and their biosynthesis pathway was also captured (Fig 2). Consequently, the general composition of enriched pathways indicates that targeted proteins serve in the pathways that are crucial for amino acid and cell membrane metabolisms.

Fig 2

Pathway enrichment analysis result.

The values on the x axes indicates number of drug targets in the related pathway with no homology to human proteins.

Druggability analysis was performed to identify potential drug targets that can be targeted with chemical molecules. The druggability analysis aims to identify proteins that have a high affinity to bind to drug-like molecules since some proteins do not have this property [18]. Here, 43 potential drug targets that have high affinity to bind drug-like molecules were determined among 89 non-homologous pathogen proteins (S3 Table). Determining potential drug targets that are broadly distributed among other harmful bacteria is one of the important steps of the prioritization process. [18]. To identify such potential drug targets, broad-spectrum analysis was performed. Finally, 28 potential drug targets that are non-homologous to human proteins, druggable and broadly distributed among other bacteria were determined. The list of final potential drug targets is given in Table 3 and S4 Table. Of 28 predicted potential drug targets, 20 were reported as essential in the DEG database. When we ranked the drug target list in terms of broad spectrum score, i.e. number of pathogenic bacteria with significantly similar sequence of the gene, 16 of the top 20 genes were essential based on DEG (80%).

Table 3

Model-derived potential drug targets that obey three criteria: No homology to human proteins, druggable, and broad-spectrum behaviour.

Locus Names	Gene Symbol	Protein Name	Pathway	Reported at DEG
STM1824	pabB	Aminodeoxychorismate synthase component 1	Folate biosynthesis	Yes
STM0087	folA	Dihydrofolate reductase	Folate biosynthesis	Yes
STM0183	folK	2-amino-4-hydroxy-6-hydroxymethyldihydropteridine pyrophosphokinase	Folate biosynthesis	No
STM3295	folP	Dihydropteroate synthase	Folate biosynthesis	No
STM0064	dapB	4-hydroxy-tetrahydrodipicolinate reductase	Biosynthesis of amino acids, L-lysine biosynthesis via DAP pathway	Yes
STM0213	dapD	2,3,4,5-tetrahydropyridine-2,6-dicarboxylate N-succinyltransferase	Biosynthesis of amino acids, L-lysine biosynthesis via DAP pathway	Yes
STM0207	mtnN	5’-methylthioadenosine/S-adenosylhomocysteine nucleosidase	Biosynthesis of amino acids, Cysteine and methionine metabolism.	No
STM3486	aroB	3-dehydroquinate synthase	Biosynthesis of amino acids, Phenylalanine, tyrosine and tryptophan biosynthesis	No
STM2384	aroC	Chorismate synthase	Biosynthesis of amino acids, Phenylalanine, tyrosine and tryptophan biosynthesis	No
STM3862	glmU	Bifunctional protein GlmU	UDP-N-acetyl-alpha-D-glucosamine biosynthesis	Yes
STM2094	rmlC	dTDP-4-dehydrorhamnose 3,5-epimerase	Polyketide sugar unit biosynthesis, Streptomycin biosynthesis	No
STM1772	kdsA	2-dehydro-3-deoxyphosphooctonate aldolase	Lipopolysaccharide biosynthesis	Yes
STM3316	kdsC	3-deoxy-D-manno-octulosonate 8-phosphate phosphatase KdsC	Lipopolysaccharide biosynthesis	No
STM0988	kdsB	3-deoxy-manno-octulosonate cytidylyltransferase	Lipopolysaccharide biosynthesis	Yes
STM0310	gmhA	Phosphoheptose isomerase	Lipopolysaccharide biosynthesis	No
STM0228	lpxA	Acyl-[acyl-carrier-protein]—UDP-N-acetylglucosamine O-acyltransferase	Lipopolysaccharide biosynthesis	Yes
STM0134	LpxC	UDP-3-O-acyl-N-acetylglucosamine deacetylase	Lipopolysaccharide biosynthesis	Yes
STM1200	tmk	Thymidylate kinase	Pyrimidine metabolism	Yes
STM1707	pyrF	Orotidine 5’-phosphate decarboxylase	Pyrimidine metabolism	Yes
STM1426	ribE	Riboflavin synthase, alpha chain	Riboflavin metabolism	Yes
STM0417	ribH	6,7-dimethyl-8-ribityllumazine synthase	Riboflavin metabolism	Yes
STM0045	ribF	Riboflavin biosynthesis protein	Riboflavin metabolism	Yes
STM3307	murA	UDP-N-acetylglucosamine 1-carboxyvinyl transferase	Peptidoglycan biosynthesis. Amino sugar and nucleotide sugar metabolism	Yes
STM0129	murC	UDP-N-acetylmuramate—L-alanine ligase	Peptidoglycan biosynthesis. D-Glutamine and D-glutamate metabolism	Yes
STM0123	murE	UDP-N-acetylmuramoyl-L-alanyl-D-glutamate—2,6-diaminopimelate ligase	Peptidoglycan biosynthesis.Lysine biosynthesis	Yes
STM0128	murG	UDP-N-acetylglucosamine—N-acetylmuramyl-(pentapeptide) pyrophosphoryl-undecaprenol N-acetylglucosamine transferase	Peptidoglycan biosynthesis	Yes
STM0124	murF	UDP-N-acetylmuramoyl-tripeptide—D-alanyl-D-alanine ligase	Peptidoglycan biosynthesis.Lysine biosynthesis	Yes
STM3725	coaD	Phosphopantetheine adenylyltransferase	Pantothenate and CoA biosynthesis	Yes

Analysis of the prioritized drug targets

The prioritized drug targets (Table 3) were clustered into pathways that are crucial for the survival of the pathogen. There are four proteins, pabB, folA, folK, folP, that take part in folate biosynthesis pathways. For nearly all organisms, the folate biosynthesis pathway is one of the indispensable pathways in order to maintain life. Folates are necessary for the production of essential biomolecules such as nucleic acids and amino acids. Most bacteria, fungi and plants can synthesize folate, while animal cells take it up from external sources [38]. This pathway is a potentially promising drug target since human cells do not have a folate synthesis mechanism that might be manipulated by a pathogen. Five of the prioritized drug targets have functionality in the biosynthesis of amino acids (dapB, dapD, mtnN, aroB, aroC). dapB and dapD take place in L-lysine biosynthesis via diaminopimelic acid (DAP) pathway, and the side product of this pathway is m-DAP, which is an essential biomolecule for peptidoglycan cell wall for gram-negative bacteria [39]. glmU has a role in the production of UDP-N-acetyl-alpha-D-glucosamine, which is essential for bacterial cell wall [40]. rmlC has an important role in the synthesis of L-rhamnose, which is an important saccharide for the virulence of some pathogens including S. Typhimurium. The absence of L-rhamnose biosynthesis pathway in human cells makes this drug target more appealing [41]. The survival of bacterium depends on the integrity of cell envelope. kdsA, kdsC, kdsB, gmhA, lpxA and LpxC have roles in the production of lipopolysaccharides, which is critical for the formation of cell envelope [42]. tmk and pyrF are involved in pyrimidine metabolism, which is crucial for all living organisms. tmk catalyzes the phosphorylation of thymidine 5’-monophosphate, which is an essential reaction for pyrimidine synthesis [43]. ribE, ribF and ribH are required for the production of riboflavin, which is a precursor of flavin mononucleotide (FMN) and flavin adenin dinucleotide (FAD). Riboflavin synthesis is a pathway required for the survival of gram-negative bacteria in the absence of external riboflavin synthesis [18]. murA, murC, murE, murG and murF are involved in the synthesis of peptidoglycans, which is an essential ingredient for bacterial cell wall biogenesis [44]. murA catalyzes the first step of peptidoglycan biosynthesis, and deletion of murA leads to death of Escherichia coli and Streptococcus pneumoniae. Fosfomycin is an antibiotic that targets murA in order to kill bacteria [45]. Mur ligases are known to be attractive drug targets because of their role in bacterial cell wall formation [46]. Pantothenate is a main precursor of coenzyme A, and its absence leads to deficiency in bacterial growth. coaD is involved in the fourth step in the coenzyme A biosynthesis pathway, which was investigated before as a suitable antibiotic target [47].

Identification of potential drugs for pabB

4-aminobenzoic acid (PABA) synthesis is an attractive antibiotic target since it is an essential biomolecule for many pathogens and it does not have a human counterpart. PABA has two main functionalities in the bacteria; (i) it is a substrate for folic acid pathway, which is critical for survival of pathogen, (ii) it is a precursor in coenzyme Q biosynthesis, which is essential for virulence [48, 49]. PABA is synthesized in two steps, and the first step is catalyzed by pabA and pabB enzymes by converting chorismate to 4-amino-4-deoxychorismate [50]. And, the second step is the production of PABA from 4-amino-4-deoxychorismate. pabB was detected as one of the putative drug targets in this study by our drug target prioritization pipeline. We specifically focused on pabB in the rest of our study as a drug target to eliminate Salmonella infections since it has critical functionality in PABA synthesis pathway, which is not represented in human cells. pabB was reported to be essential experimentally in the DEG database, and, among the identified drug targets with experimental validation (Table 3), it ranks 7th in terms of the number of pathogenic bacteria strains that carry a gene with high sequence similarity based on our broad-spectrum analysis. Additionally, we investigated the importance of pabB in the pathogen-host integrated GMN model. Interactions of the metabolites of the pabB reaction (chorismate, 4-amino-4-deoxychorismate, L-glutamate, L-glutamine) with the metabolites of other reactions were visualized by creating a metabolite-metabolite interaction network (S1 Fig). L-glutamate is directly related to numerous amino acids such as alanine, leucine and asparagine. On the other hand, 4-amino-4-deoxychorismate is indirectly related with the production of thymidine through tetrahydrofolate. It is also linked with the production of adenine through R-Pantoate. Therefore, DNA synthesis is dependent on the production of 4-amino-4-deoxychorismate through adenine and thymine synthesis which cannot be synthesized when pabB is inhibited.

There is a very similar protein to this putative target in Escherichia coli, which is also called pabB. Formic acid, which is widely used as an antibacterial agent in fodders, was reported in DrugBank as a compound that targets pabB in Escherichia coli (strain K12) [51]. To our knowledge, pabB was not offered or examined as a drug target for Salmonella species before. Hereby, protein docking and molecular dynamic analysis were performed to determine novel molecules that can inhibit S. Typhimurium growth by binding pabB. The protein Aminodeoxychorismate synthase component 1, which belonged to S. Typhimurium (strain LT2 / SGSC1412 / ATCC 700720) with a UniProt ID: P12680 (PABB_SALTY), was taken for elaborative structural and functional studies. The 3D structure of the protein is the core requirement to carry out protein molecular docking studies and its unavailability necessitated the modeling of the three-dimensional structure. Herein, first step was the comparative 3D modeling of the query protein with the application of MODELLER Software. The crystal structure of 4-amino4-deoxychorismate (ADC) synthase (PDB ID: 1K0E) was chosen as a template as its percent identity was 76.38% and query coverage was 99% against P12680. The MODELLER yielded five 3D structures of P12680 protein, out of which the best model with higher accuracy was selected (S2 Fig) after thorough quality assessment using PROCHECK server. The Verify3D passed the modelled structure with 87% indicating 80% amino acids of the built model having score > = 0.2 in the 3D/1D profile. Further, the Ramachandran plot (S3 Fig) showed 92.7% residues in the most favoured region whereas only 0.3% residues were reported in disallowed region. Once the model was finalized, Mg2+ ion was incorporated into the built 3D structure as magnesium ion is reported as the cofactor within the target protein in the UniProt indicating its role in the catalysis activity. Therefore, Mg2+ ion was complexed near the active site residues reported in UniProt and in literature [48].

The DoGSiteScorer produced 10 probable binding pockets, out of which pocket number one was selected as the binding site to utilize for molecular docking. The selected pocket has a druggability score of 0.83, suggesting the prime region for drug binding. The binding pocket was analyzed in Chimera software [52], and validated as it covers the active site amino acid residues reported in the UniProt along with the amino acid residues coordinating with the Mg2+ ion. The predicted binding site can be seen in (S4 Fig).

Once all the 54,000 drug-like compounds were screened against the target protein (P12680), the lowest binding energy conformation of each single 54,000 compound(s) was obtained. The overall binding energies ranked against the P12680 are represented through a histogram in Fig 3.

Fig 3

Histogram illustration of the overall binding energy retrieved through virtual screening.

The thorough analysis of the results suggested that 1659 compounds showed promising binding free energy ranging from -12.42 kcal/mol to -9.04 kcal/mol, while 22,231 compounds having binding free energies within the range of -5.66 kcal/mol to -2.28 kcal/mol. Moreover, top ten compounds that were docked near the active site of the target protein were (having binding free energy of -12.42, -12.02, -11.92, -11.90, -11.91, -11.98, -11.77, -11.73, -11.42, -11.72) retrieved to further investigate their chemical interactions with amino residues of target protein. The aforementioned top ten compounds with their ZINC IDs are reported in (S5 Table).

The protein-ligand complex of each top ten compound was explored to inspect the chemical bonds and interactions occurring within the protein-ligand complex. LigPlot+ generates the atomic interactions taking place among the target protein and ligand (drug) through hydrogen bonds and hydrophobic contacts. Each of the ten protein-ligand complex analyzed in LigPlot+ has shown ligand interactions with the Mg2+. The reported active site residues from UniProt (i.e. Lys275 and Glu259) were found to be interacting with the top ten ligands, out of which Lys275 can be seen making hydrogen bond with five ligands and hydrophobic contacts with two ligands. (Table 4). On the other hand, Glu259 makes hydrophobic contacts with six of the ligands. The most common amino acid residues interacting with top ten ligands through hydrogen bonding were identified as Thr277, Gly427, Glu440, Lys444, Lys275, Arg411 and Glu259. Finally, each of the amino acid residues making hydrogen bonds and hydrophobic contacts interaction with the top ranked 10 compounds is presented in Table 4 and visualized in S5–S7 Figs. There is not any report in literature about the use of these compounds against pathogenic bacteria. Therefore, they remain open to experimental validation.

Table 4

Protein residues involved in the hydrogen and hydrophobic interactions with top ten best ranking compounds (ligands), analyzed through LigPlot+.

S.No	ZINC IDs	Residues making hydrogen bond interaction	Residues making Hydrophobic contacts
1	ZINC7879733	Thr277(A), Gly427(A) Lys444(A), Arg411(A) Lys275(A)	Ile410(A), Ala424(A), Gly427(A) Asn214(A) Trp391(A), Ser423(A), Gly426(A), Glu259(A) Ile368(A), Cys422(A), Gly425(A), Gly276(A), Ile274(A)
2	ZINC15179659	Gly427(A), Arg411(A), Lys444(A), Thr277(A), Glu440(A)	Lys275(A), Trp391(A), Val445(A), OThr277(A) Ala424(A), Ser367(A), Ile410(A), Gly425(A) Ser423(A), Ile368(A), Ile448(A), Asn214(A), Gly426(A) Thr412(A), Cys422(A), GIle274(A), ly276(A)
3	ZINC14880941	Thr277(A), Ser423(A) Gly427(A), Gly276(A) Glu440(A), Trp391(A) Asn214(A)	Lys444(A), Arg411(A), Val445(A), Ile410(A) Lys275(A), Ile274(A), Gly425(A), Ala424(A), Ile368(A) Gly426(A), Ile448(A), Glu259(A), Cys422(A)
4	ZINC58542694	Arg411(A), Thr277(A), Glu440(A), Lys275(A), Gly427(A),	Lys444(A), Gly276(A), Trp391(A), Ile410(A), Ser423(A), Ala424(A), Ile274(A), Asn214(A), Ile368(A), Val445(A), Gly425(A), Gly426(A), His340(A), Ser367(A),
5	ZINC1201089024	Trp391(A), Lys275(A), Asn214(A), Arg411(A), Glu440(A), Thr277(A), Gly427(A),	Lys444(A), Gly276(A), Ile410(A), Ile274(A), Gly425(A), Ile368(A), Val445(A), Glu259(A), Ser423(A), Gly426(A) Cys422(A), Ile448(A),
6	ZINC27071723	Lys444(A), Trp391(A), Lys275(A), Glu440(A), Thr277(A), Gly427(A),	Arg411(A), Gly276(A), Ile274(A), Val445(A), Ser423(A), Gly425(A), Ile410(A), Cys422(A), Gly426(A), Thr412(A), Ala424(A), Ile448(A), Asn214(A), Glu259(A), Ile368(A),
7	ZINC7133393	Arg411(A), Gly425(A), Glu440(A), Thr277(A), Gly427(A)	Lys444(A), Gly276(A), Trp391(A), Ile410(A), Ala424(A), Val445(A), Ser423(A), Glu259(A), Asn214(A), Gly426(A)Ser367(A)
8	ZINC7879735	Gly427(A), Lys444(A), Arg411(A), Glu440(A), Thr277(A),	Lys275(A), Gly425(A), Gly276(A), Asn214(A), Ile274(A), Ala424(A), Ile410(A), Ser367(A), Trp391(A), Ser423(A), Gly426(A), Cys422(A), Val445(A),
9	ZINC58542238	Gly427(A), Gly425(A), Glu440(A), Lys275(A), Thr277(A)	Arg411(A), Lys444(A), Gly276(A), Trp391(A), Ser423(A), Val445(A), Ile274(A), Ile410(A), Asn214(A), Ala424(A), Gly426(A), His340(A), Ile368(A), Thr412(A), Cys422(A), Ile448(A),
10	ZINC7538530	Arg411(A), Glu440(A), Thr277(A), Gly427(A),	Lys444(A), Gly276(A), Lys275(A), Ile274(A), Gly425(A), Trp391(A), Ile410(A), Ser423(A), Ala424(A), Gly426(A), Asn214(A), Glu259(A), Cys422(A), Val445(A)

Conclusions

Analysis of pathogen-host integrated GMN models is a well-suited approach in terms of identifying novel drug targets by considering pathogen-host interactions. It allows tracking the response of pathogen and host simultaneously during infection by mapping infection-induced dual-transcriptome data. There are some pathogen-host GMN models published [23, 24], but to our knowledge, this is the first study in the literature that determines drug targets by analysing condition-specific pathogen-host integrated GMNs created by mapping dual-transcriptome data. We here reconstructed and analyzed condition-specific integrated GMN models to identify novel drug targets for S. Typhimurium induced infections. We used prioritization steps to identify best suitable novel drug targets. After prioritization processes, we identified 28 putative drug targets, and pabB was chosen as a high ranked drug target based on our prioritization pipeline, literature information and novelty. Subsequently, homology and molecular docking analyses were performed to identify candidate compounds that inhibit pabB. The top ten compounds in terms of binding free energy were identified and reported. Consequently, analysing S. enterica metabolism inside the host cell has enabled us to comprehend metabolic alterations in both HeLa cells and S. enterica along with determining novel drug targets. Future studies may provide more arguments for the proposed drug targets and inhibitors in this study. This study can be used as a guideline for creating and analysing condition specific pathogen-host GMN models.

Materials and methods

The flowchart of the pipeline followed in this study is given in Fig 4. Each step is detailed in the sections below.

Fig 4

The flowchart of the pipeline followed in this study.

Transcriptome data

The dual RNA-seq data of infected HeLa cells and S. Typhimurium strain SL1344 [25] was downloaded from NCBI Gene Expression Omnibus (GEO) Database [53]. The dataset ID in the GEO database is GSE117236. The samples collected at the beginning of infection and the post-infection data at 8th and 16th hours were used in this study. Each time point included duplicate samples. Principal component analysis (PCA) was used to identify any possible outliers in the data, and no outliers were detected (S8 Fig).

Pathogen-host integrated genome scale metabolic network

Two different genome-scale metabolic models from the literature were used to reconstruct an integrated pathogen-host genome scale metabolic network (GMN). A genome-scale metabolic network of S. Typhimurium, called stm_v1.0, consisting of 2,545 reactions controlled by 1,271 genes was used as the pathogen metabolic network [21]. As the host GMN, a recent reconstruction of human metabolism with a substantial amount of curation, called iHsa, was used, which covered 8,336 reactions and 2,315 genes [54]. S. Typhimurium is an intracellular pathogen. Therefore, the GMN of the pathogen was placed inside the cytoplasm of the host GMN as a separate compartment to create the pathogen-host integrated GMN. Here, an extensive literature research was performed to identify cytosolic host metabolites that can be consumed by the pathogen during infection. An exchange reaction with the extracellular environment must be available for these metabolites in the S. Typhimurium metabolic network. 38 such metabolites were identified [26, 29, 55–57], and they were allowed to be taken up by the S. Typhimurium GMN from the cytoplasm of the host (S6 Table). In addition, the pathogen was allowed to secrete all its exchange metabolites to the host cytoplasm as defined by the model secretion reactions (S7 Table). As a result, the pathogen-host integrated GMN consisting of 11,029 reactions controlled by 3,586 genes was created by using COBRA Toolbox v.3.0 on MATLAB programming platform [58].

Flux balance analysis

Pathogen-host integrated GMN was analyzed using flux balance analysis (FBA) method to predict fluxes associated with the infection times studied. FBA searches the solution space defined by mass balance and reaction reversibility constraints to find an optimal solution with the help of an objective function. FBA assumes that the system is at steady state, i.e. the concentrations of intracellular metabolites do not change over sufficiently long time, leading to linear mass balance constraints [59]. The objective function of pathogen-host integrated GMN was defined based on a weighted relationship between host and pathogen biomass, and the biomass composition formulas were taken from iHsa and stm_v1.0 models (Eq 1) [21, 54, 60]. HB, PB and TB are host biomass, pathogen biomass and total biomass respectively in Eq (1), where α and β are maximum host and pathogen biomass production rates in order. Using the constructed pathogen-host GMN, α and β were calculated first with FBA via maximization of HB and PB reactions separately. Later, maximum host and pathogen biomass production rates (α and β), which were calculated as 0.39 and 0.28 respectively, were added to the equation as weights of HB and PB to get TB, which represents the balanced effect of HB and PB. Eq 1 was added as a reaction to the pathogen-host integrated GMN and later set as the objective function.

The upper bound of glucose uptake rate of S. Typhimurium was set to 5 mmol/gDW/h in simulations based on the studies of Thiele and her coworkers [21]. The maximum uptake rate of other available carbon sources for S. Typhimurium inside the host cell was limited to 20% of its glucose uptake rate. Oxygen uptake rate of S. Typhimurium was set to 1 mmol/gDW/h to mimic the hypoxic environment during infection [61]. The host cell was allowed to utilize only metabolites that were found in the Dulbecco’s Modified Eagle’s Medium (DMEM) since the infection experiment for the dual RNA-seq data, which was used in creating condition-specific integrated GMN models, was carried out in this medium (S8 Table). In GMN models, alternate optima can be an issue in interpreting the results of FBA since there might be multiple flux distributions that result in the same value for the objective function. Minimization of the sum of squares of all flux values was applied to prevent alternate optima [62]. The principle of this method was proposed based on the accomplishment of the cellular goals with minimal resource expenditure since the flux values are an indication of the amount of depletion of resources [63].

Integrating transcriptome data with pathogen-host integrated GMN

Condition specific GMNs were generated by mapping dual RNA-seq data on the pathogen-host GMN to simulate infection states at different time points. Gene Inactivity Moderated by Metabolism and Expression (GIMME) algorithm was used as the mapping algorithm to generate condition specific GMNs [22]. GIMME determines active reactions based on a threshold put on the mRNA levels in the data, where reactions below the thresholds are set as inactive. GIMME generates a GMN with the desired functionality using the objective fraction parameter, and it adds the reactions in the inactive set back if their removal affects the desired functionality [22]. The threshold value was determined as the quarter of the mean of the transcriptome data. Since the average gene expression values were much higher in HeLa cell compared to S. Typhimurium in the utilized dual RNA-seq data (S9 Fig), the threshold value was separately determined for both organisms. Then, since GIMME algorithm accepts a single threshold, the difference between the organism-specific threshold values were added to the reaction scores of S. Typhimurium, and the threshold value obtained from the human transcriptome was used as the threshold value in GIMME simulations. The objective fraction parameter was set to 0.2 in order to ensure that the condition specific integrated GMN produces at least 20% of the maximum TB. Three different condition specific GMNs were produced as a result to represent the start of infection (0th hour) and infections at 8th and 16th hours.

Identification of drug targets

Gene deletion analysis is a widely used approach in constraint-based metabolic modeling to predict potential drug targets, and it is performed by in silico deletion of genes in the GMN [10]. The analysis aims to obtain essential genes for the desired functionality of a GMN, such as preventing growth of pathogen. FBA can be used to predict essential genes in an organism by setting the rate of the associated reaction(s) to zero for each gene. If inactivation of the reaction(s) lead to zero growth rate, the gene is essential for the pathogenic organism and it can be used as a drug target. In this study, gene deletion analysis was performed to identify potential drug targets that can restore the metabolic changes in the host cell caused by S. Typhimurium induced infection. GIMME-based condition-specific GMNs were used in the analysis. A nonzero rate for the production of HB together with zero rate for PB was used as the desired functionality of the condition specific GMNs in gene deletion analysis. Using the FBA technique to maximize HB and PB separately, essential genes were obtained, and the enzymes that catalyze these reactions were chosen as potential drug targets.

The identified drug targets were further filtered based on the approach applied elsewhere [18]. Briefly, the steps in the approach followed are as follows: (i) Homology analysis is a critical part of drug target selection process, and the aim of the analysis is determining drug targets that are not similar to host proteins in order to avoid side effects of drugs. Homology analysis was performed using BLASTp algorithm, and drug targets that are not similar to human proteins were identified [64]. As BLASTp parameters, cut off for the expected value (E-value) was chosen as 1x10-4, and the maximum sequence identity for the determination of human-non-homolog drug targets were chosen as 30% [18]. Pathway enrichment analysis was performed with the identified human-non-homolog drug targets by using KOBAS V3.0 [65] to elucidate pathways that are expected to be crucial for the survival of S. Typhimurium inside the host. (ii) Another important step in the model-based drug target selection process is the determination of druggable proteins, and the goal of this selection is identifying drug targets that can be targeted with drug-like chemical compounds. Druggable targets were identified using BLAST algorithm in Drugbank database [51], and the E-value was chosen as 10−25 [18]. (iii) Broad spectrum analysis was performed for the identification of drug targets that are broadly distributed among other bacteria. Broad spectrum analysis is very beneficial in order to determine drug targets that can be effective against co-infections or multiple infections. In addition, the proteins that are broadly distributed among other bacteria may indicate low mutation rate, so developing antibacterial resistance can be harder for the targeted bacteria. Broad spectrum analysis was performed using PBIT web browser [66], which contains protein sequences of 181 pathogenic organisms. E-value cut-off of 1x10-5, bit score of 100 and sequence identity of 35% were chosen as parameter values [18, 66]. For the broad spectrum target determination criteria, targets available in at least 40 pathogenic strains were chosen [18].

Homology modeling of target protein and active / binding pocket prediction

As the 3D structure of target protein was not available in the Protein Data Bank (PDB) [67], the homology modeling was performed to model the target protein’s 3D structure. In the first step, the appropriate template protein was searched by BLASTp against the PDB database. The template that matches the criteria of query coverage ≥ 90% and percent identity ≥ 70% was selected. The Modeller software v9.20 was used to perform comparative protein structure modelling [68-71]. The Modeller works by satisfying the spatial restraints along with employing specific geometrical calculations generating possible coordinate(s) for the location of each single atom of target protein [70]. The homology modeling was performed through Modeller by implementing a series of python script(s), which resulted in generating five models with their Discrete Optimized Potential Energy (DOPE) values. The 3D model with the lowest DOPE value was selected as the predicted 3D structure of the target protein. The predicted 3D model was then validated through verify3D and Ramachandran Plot via PROCHECK server to ensure its reliability and quality. The Verify3D program was used to interpret the quality of the built model, as Verify3D computes the compatibility of 3D model of target protein against its amino acid sequence [72]. Further, Ramachandran Plot was analyzed to assess the stereo-chemical properties of the predicted 3D model [73].

The built target protein model was subjected to DoGSiteScorer to identify potential binding pockets. The DoGSiteScorer is an automated grid-based program which utilizes difference of Gaussian filter to identify the potential binding pocket present within the protein [74, 75]. The program provides ten potential binding pockets with druggability scores. Out of the ten predicted pockets, the one with the highest druggability score is supposed to be a rich binding pocket and therefore, selected for study. The program also evaluates the depth, surface area and volume for each predicted binding pocket.

Docking-based virtual screening

To carry out the rigorous virtual screening of drug-like molecules against the target protein, a library of chemical compounds was curated. The ZINC15 database was used to retrieve the drug-like molecules following certain criteria of compound’s molecular weight and logP values [76]. The drug-like compounds having logP value ≤ 5 and Molecular Weight ≤ 375 Daltons were obtained from the database in SDF format. A total of 54,000 drug-like compounds were compiled and prepared into the required PDBQT format by using Open Babel. All the 54,000 compounds were minimized with force field MMFF94 via Open Babel [77].

The AutoDock-GPU [78] was chosen to execute molecular docking and virtual screening of curated drug-like molecules (ligands) library against the target protein (P12680). The AutoDock GPU was chosen because of the prolonged execution times whilst using AutoDock4. AutoDock-GPU, which is an OpenCL and cuda based implementation of Autodock4, was created to utilize large number of GPU cores and speed up docking by using parallel processing [79, 80]. So, to execute molecular docking, the target protein’s modelled 3D structure was prepared by adding hydrogens and its conversion into the required format of PDBQT. Afterwards, the grid search box dimensions were set carefully to cover the predicted binding site retrieved from the DoGSiteScorer. Further docking steps were carried out to create docking parameter files, and, eventually, the grid maps FLD file, docking parameters files, GPF and DPF were made. Once the necessary files were created, the molecular docking was performed to screen all the 54,000 compounds with 20 Genetic Algorithm (GA) runs against target receptor’s binding site. The virtual screening yielded binding free energy for 20 runs of all the 54,000 compounds. Ultimately, the lowest binding energy of each compound was extracted for structural and functional evaluation. Moreover, re-docking step was performed by execution of all the Autodock4 steps utilizing only one ligand against the target protein to validate predicted binding site.

LigPlot+ program was used to generate 2D diagrams of the target protein and ligand complex. For the LigPlot+ analysis, protein-ligand complex in PDB format was taken. The atomic interaction(s) within the diagram is shown, where ligand and protein’s interactive residues are represented in ball-stick format. From the diagram, the amino acid residues of the target protein making chemical interactions with the ligand can be identified. The LigPlot+ highlights the hydrogen bonding within the atoms of protein and ligand with green dotted lines.

Metabolic network of pabB related metabolites.

The red arrows indicate the reaction controlled by pabB (L-Glutamine + chorismate -> L-Glutamate + 4-amino-4-deoxychorismate).

(TIF)

Click here for additional data file.

The predicted 3D structure of P12680.

(TIF)

Click here for additional data file.

Ramachandran plot of the predicted 3D structure of target protein (P12680).

(TIF)

Click here for additional data file.

The predicted binding pocket is shown as surface in sand brown color while the protein chain as ribbon in blue color.

From the figure, it can be seen that Mg2+ ion is submerged within the predicted binding pocket.

(TIF)

Click here for additional data file.

A) is representing the LigPlot+ of P12680 showing atomic interaction between protein (residues/ Mg2+ ion) and ligand (ZINC7879733). The atomic linkages due to hydrogen bonding can be identified from the diagram. Similarly, B) represents the LigPlot+ analysis of P12680 and ligand (ZINC15179659), C) represents the LigPlot+ for P12680 and ligand (ZINC14880941), and D) represents the LigPlot+ for P12680 and ligand (ZINC58542694).

(TIF)

Click here for additional data file.

A) is the LigPlot+ for P12680 and ligand (ZINC1201089024), all the atomic linkages occurring between the protein-ligand complex can be analyzed. Likewise, B) represents the LigPlot+ for P12660 and ligand (ZINC27071723), C) is LigPlot+ for P12680 and ligand (ZINC7133393) and D) is LigPlot+ for P12680 and ligand (ZINC7879735).

(TIF)

Click here for additional data file.

A) is the LigPlot+ for P12680 and ligand (ZINC58542238) complex and B) is the LigPlot+ for P12680 and ligand (ZINC7538530) complex. The overall amino acid residues of P12680 which interact with each of the top ten compounds (ligands) making hydrogen bonds and hydrophobic contacts.

(TIF)

Click here for additional data file.

Principal component analysis (PCA) of GSE117236.

The red, blue and green dots represent beginning of infection, 8th hour of infection and 16th hour of infection respectively.

(TIF)

Click here for additional data file.

Box plots of gene expression values of HeLa cell and S. Typhimurium.

(TIFF)

Click here for additional data file.

140 essential genes based on DEG database.

(DOCX)

Click here for additional data file.

89 potential drug targets that are not similar to human proteins.

(DOCX)

Click here for additional data file.

Potential drug targets that have high affinity to bind drug-like molecules.

(DOCX)

Click here for additional data file.

The list of final potential drug targets.

The last column indicates broad spectrum analysis results.

(DOCX)

Click here for additional data file.

Binding energy, Zinc IDs, 1D and 2D structure of each of top 10 compounds.

(DOCX)

Click here for additional data file.

Metabolites that can be consumed by S. Typhimurium inside host cytoplasm in pathogen-host GMN model.

(DOCX)

Click here for additional data file.

The exchange metabolites of S. Typhimurium that match with the cytoplasmic metabolites of human cell.

(DOCX)

Click here for additional data file.

DMEM medium constraints for the host GMN used in metabolic model simulations.

(DOCX)

Click here for additional data file.

Long names of metabolites given in S1 Fig, the figüre that reports metabolite-metabolite interactions around pabB catalyzed reaction.

(DOCX)

Click here for additional data file.

The MATLAB codes to generate and analyze pathogen-host integrated genome scale metabolic models of S. Typhimurium and human.

(ZIP)

Click here for additional data file.

Transfer Alert

This paper was transferred from another journal. As a result, its full editorial history (including decision letters, peer reviews and author responses) may not be present.

17 Feb 2022

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Roman G. Gerlach

Academic Editor

PLOS ONE

Journal Requirements:

When submitting your revision, we need you to address these additional requirements.

1. Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Thank you for stating the following in the Acknowledgments Section of your manuscript:

“This work was supported by TUBITAK, The Scientific and Technological Research Council of Turkey (Project Code: 316S005) and by PSF, The Pakistan Science Foundation [Project Code: PSF-TUBITAK/S-HEJ (04)].”

We note that you have provided additional information within the Acknowledgements Section that is not currently declared in your Funding Statement. Please note that funding information should not appear in the Acknowledgments section or other areas of your manuscript. We will only publish funding information present in the Funding Statement section of the online submission form.

Please remove any funding-related text from the manuscript and let us know how you would like to update your Funding Statement. Currently, your Funding Statement reads as follows:

“This work was supported by TUBITAK, The Scientific and Technological Research Council of Turkey (Project Code: 316S005) and by PSF, The Pakistan Science Foundation [Project Code: PSF-TUBITAK/S-HEJ (04)]. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.”

Please include your amended statements within your cover letter; we will change the online submission form on your behalf.

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: I Don't Know

Reviewer #2: I Don't Know

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: In this manuscript, Kocabas et al took advantage of the available dual RNA-seq data from Salmonella-infected cells, and constructed in silico an integrated metabolic network for Salmonella and host cells. Based on the metabolic network, the authors predicted 140 essential bacterial metabolic genes and further identified pabB as a potential druggable target. In silico docking analysis suggested PabB might interact with small molecules for drug development.

Though most of the analysis were purely in silico with little credibility, this is an innovative approach and an interesting idea in my personal opinion. The study provides additional angles to decode other RNA-seq data and may help promote the understanding of Salmonella infection processes.

Reviewer #2: Kocabaş et al report the reconstruction of transcriptome-specific pathogen-host integrated models using dual RNA-seq data of Salmonella Typhimurium and infected HeLa cells, which predicts pabB gene of the pathogen as one of the top druggable targets. The authors have also performed structure-based virtual screening of 54,000 compounds against pabB to identify 1659 compounds, leading to a final top 10 best-ranking compounds. The manuscript highlights the use of dual RNA seq based metabolic modeling as an approach to identify drug targets - which has scientific merit - however it would benefit by considering some experimental validations and additional performance assessment of the model predictions. The host (HeLa cell) metabolic model, iHsa (Blais et al, 2017) and the pathogen model, stm_v1.0 (Thiele et al 2011) are already published ones and the concept of pathogen-host integration is also reported already for Mycobacterium (Bordbar et al 2010, Rienksma et al 2019). The novelty of this work is the use of dual RNA seq to constrain the metabolic model, that is based on GIMME algorithm (Becker et al 2008), which is well-known for developing cell-specific models. Considering the overlap between reported concepts, it requires rigorous performance assessments to make it condition-specific and to be considered as the first application to Salmonella.

Comments:

• The analyses of the model performance could have been more extensive while comparing with experimental validations (For e.g., only secretion of 3 metabolites are compared in the manuscript – more such predictions need to be validated - in the host-side of the network as well).

• The model assessment for precision and recall, to deduce its accuracy of gene essentiality prediction needs to be performed in comparison with the experimental datasets (e.g., mutants in Salmonella for the identified gene sets?) or literature evidence to categorize the predictions as “true positives” and “false positives”. Among the 140 gene targets from gene essentiality predictions, how many of them have literature evidence as essential for Salmonella inside host (HeLa) cells or are all novel predicitons?

• The prioritization criteria for selecting pabB from 28 potential drug targets - that are non-homologous to human proteins, druggable and broadly distributed among other bacteria, is only based on the importance for PABA. This should be justified by mentioning how they arrived at one candidate gene based on the model predictions, for e.g., effect in biomass within host (% inhibition) or participation in higher number of essential reactions in the network?

• Metabolic modeling could be used to track the mechanism of a knockout phenotype. Could the authors sketch out the mechanism by which pabB is becoming essential in Salmonella? via metabolic network – delineating the reactions involving pabB and the effect of deleting pabB in Salmonella? How many metabolic reactions in pathogen network are linked to pabB and/or other gene targets identified? How many of these reactions become essential reactions for the pathogen within host?

• Line 383, Equation 1 - α and β are calculated via FBA and then used as weights. Will it be limited to the constraints used each time the FBA is run or is fixed for all simulations? Instead, could FBA with two objective functions (HB and PB) at the same time, by changing objective function weights would be more appropriate?

• Top 10 compounds identified and listed in Table 4, requires more explanation on the status of their current applications – any homology with natural products, or in clinical trial or used in other bacterial screening etc., in order to increase the application of these identified compounds.

• Experimental screening for at least 3 compounds is required to support the claim that it can be a useful drug against Salmonella inside host cell. For e.g., Treating the Hela cells infected with Salmonella in vitro should inhibit the growth of the pathogen as predicted. Or finding the binding specificity of these compounds for pabB gene.

• Figure 1 would benefit by having two panels – separately for A) host and B) pathogen and having the column bar side by side to represent 0th and 16th hour instead of superimposing them.

• The dataset used also has an 8h time point. Why was it excluded from the analysis?

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

4 Apr 2022

We have updated the manuscript based on the editorial comments by removing Acknowledgement section and by making it compatible with PLOS One formatting styles. We include our responses to reviewer comments below.

Replies for Reviewer Comments

Reviewer 2:

1) The analyses of the model performance could have been more extensive while comparing with experimental validations (For e.g., only secretion of 3 metabolites are compared in the manuscript – more such predictions need to be validated - in the host-side of the network as well).

We agree with the reviewer that additional tests on the model performance would strengthen the manuscript. On the other hand, the studies that report measured reaction fluxes of host and pathogen during infection of S. Typhimurium are very limited in the literature. We could only find information about the secretion rates of these three pathogen metabolites to evaluate the performance of the model. To provide extra validation of the model -as recommended by the reviewer in the next comment- we additionally compared our gene essentiality predictions with the literature (details given in the next reply). We believe the metabolite secretion prediction and gene essentiality prediction performances provide evidences on the suitability of our pathogen-host integrated genome-scale metabolic network.

2) The model assessment for precision and recall, to deduce its accuracy of gene essentiality prediction needs to be performed in comparison with the experimental datasets (e.g., mutants in Salmonella for the identified gene sets?) or literature evidence to categorize the predictions as “true positives” and “false positives”. Among the 140 gene targets from gene essentiality predictions, how many of them have literature evidence as essential for Salmonella inside host (HeLa) cells or are all novel predictions?

We are grateful to the reviewer for this comment, which enabled us to provide an additional source of validation for our pathogen-host integrated genome-scale metabolic model. Predicted essential genes were checked by using Database of Essential Genes (DEG), which reports data from three experimental gene deletion studies from rich medium experiments (Barquist et al., 2013; Khatiwara et al., 2012; Knuth et al., 2004). The results were provided in the revised manuscript under the section “Identification of Potential Drug Targets”, first paragraph and 3rd paragraph. Data were available for 137 of 140 predicted essential genes, 93 of which were reported as essential genes in at least one study (68%). Regarding the prioritized list of 28 essential genes, 20 were reported as essential in the database. When we ranked the prioritized list in terms of broad spectrum score (number of pathogenic bacteria with significantly similar sequence of the gene), 16 of the top 20 genes were essential based on the DEG database (80%). Considering that the experiments were performed in rich medium conditions, with no host cells involved, we believe the predictions provide additional validation of the pathogen-host integrated genome-scale metabolic network. We have also added a new column to Table 3, indicating whether the predicted drug target is reported to be essential in DEG.

3) The prioritization criteria for selecting pabB from 28 potential drug targets - that are non-homologous to human proteins, druggable and broadly distributed among other bacteria, is only based on the importance for PABA. This should be justified by mentioning how they arrived at one candidate gene based on the model predictions, for e.g., effect in biomass within host (% inhibition) or participation in higher number of essential reactions in the network?

All of the proposed candidate drug targets were selected based on their effect on biomass production. Inhibition of all the proposed 28 candidate drug targets leads to blocking of biomass production (zero growth rate). pabB is also reported as an essential gene in DEG database. 20 of 28 candidate targets were experimentally verified in the DEG database. When we rank those 20 genes in terms of broad spectrum behaviour, pabB is 7th in terms of number of pathogenic bacteria strains that carry a gene with high sequence similarity. The reaction associated with pabB is related to production of nucleic acids thymine and adenine as well as production of numerous amino acids in the model. Therefore, we chose pabB for docking analysis. The network that shows the relationship of pabB with critical metabolites was given as Supplementary Figure S1 in the revised manuscript (please also see our reply to the next comment). Those explanations are now included in the revised manuscript, Section “Identification of potential drugs for pabB”, as given below. We thank the reviewer for helping us clearing up the ambiguity associated with this issue in our manuscript.

“Additionally, we investigated the importance of pabB in the pathogen-host integrated GMN model. Interactions of the metabolites of the reaction (chorismate, 4-amino-4-deoxychorismate, L-glutamate, L-glutamine) with the metabolites of other reactions were visualized by creating a metabolite-metabolite interaction network (S1 Fig.). L-glutamate is directly related to numerous amino acids such as alanine, leucine and asparagine. On the other hand, 4-amino-4-deoxychorismate is indirectly related with the production of thymidine through tetrahydrofolate. It is also linked with the production of adenine through R-Pantoate. Therefore, DNA synthesis is dependent on the production of 4-amino-4-deoxychorismate through adenine and thymine synthesis, which cannot be synthesized when pabB is inhibited.”

4) Metabolic modeling could be used to track the mechanism of a knockout phenotype. Could the authors sketch out the mechanism by which pabB is becoming essential in Salmonella? via metabolic network – delineating the reactions involving pabB and the effect of deleting pabB in Salmonella? How many metabolic reactions in pathogen network are linked to pabB and/or other gene targets identified? How many of these reactions become essential reactions for the pathogen within host?

We thank again the reviewer for this comment. We investigated the effect of inhibition of pabB through metabolite-metabolite interaction network, which is indeed a feature implemented in COBRA Toolbox by one of the authors (K. Kocabas). The network is given as a supplementary figure for the revised manuscript. The results and comments were provided in the revised manuscript under the section “Identification of potential drugs for pabB”. Briefly, we tracked the metabolites that were affected from to function of pabB. Since the associated reaction is related to several amino acids as well as two nucleic acids, its blocking inhibits macromolecular synthesis, leading to zero growth.

5) Line 383, Equation 1 - α and β are calculated via FBA and then used as weights. Will it be limited to the constraints used each time the FBA is run or is fixed for all simulations? Instead, could FBA with two objective functions (HB and PB) at the same time, by changing objective function weights would be more appropriate?

We used fixed α and β in all simulations. We deliberately did not change the objective function since then it would be ambiguous to identify the source of change..

6) Top 10 compounds identified and listed in Table 4, requires more explanation on the status of their current applications – any homology with natural products, or in clinical trial or used in other bacterial screening etc., in order to increase the application of these identified compounds.

There are not any literature reports for the biological activities of the 10 identified compounds. However, the shortlisted compounds are drug-like as they obey the Lipinski rules. In addition, these compounds were derived from drug-like library i.e. the ZINC15 database. Moreover, to focus on only drug-like compounds, we further filtered the compounds in the database and only obtained the compounds having logP value ≤ 5 and Molecular Weight ≤ 375 daltons from this database, which is stated in Materials & Methods section.

7) Experimental screening for at least 3 compounds is required to support the claim that it can be a useful drug against Salmonella inside host cell. For e.g., Treating the Hela cells infected with Salmonella in vitro should inhibit the growth of the pathogen as predicted. Or finding the binding specificity of these compounds for pabB gene.

We agree that experimental validation would strengthen our theoretical findings of 10 candidate compounds. On the other hand, our research groups do not have any experimental lab facility, we are purely computational labs. We believe our major contribution to the scientific field with this work is (i) to provide a validated pathogen-host genome-scale metabolic network for S. Typhimurium, (ii) to demonstrate the integrative analysis of dual RNA-Seq data with pathogen-host genome-scale metabolic models. We have included in the revised manuscript a statement that the compounds need experimental validation (Last sentences of section “Identification of potential drugs for pabB”.

On the other hand, the following paper from our group has been well cited as well as found experimental proof within one month of this article published online:

Ref#1: Identification of chymotrypsin-like protease inhibitors of SARS-CoV-2 via integrated computational approach, Journal of Biomolecular Structure and Dynamics, (2021), 39:7, 2607-2616, DOI: 10.1080/07391102.2020.1751298

We believe our theoretical findings will still help other researchers in drug development field, and will aid in the efforts to design novel drugs.

8) Figure 1 would benefit by having two panels – separately for A) host and B) pathogen and having the column bar side by side to represent 0th and 16th hour instead of superimposing them.

We thank again the reviewer for this comment. We recreated the Figure 1 as panel A) and B). We used bidirectional bar chart instead of creating a bar chart with column bar side by side to represent 0th and 16th since it otherwise led to a too big figure that contained a lot of empty spaces.

9) The dataset used also has an 8h time point. Why was it excluded from the analysis?

We wanted to compare the beginning of the infection with the late stage of the infection. We could not see any considerable differences between 0 time point and 8h time point in terms of predicted flux values. Probably the infection does not yet show its damage to the cells in this time point. Therefore, we did not include 8h time point. We added this information in the revised manuscript under the Materials&Methods section “Transcriptome data”.

Submitted filename: response_to_reviewers.pdf

Click here for additional data file.

4 May 2022

5 May 2022

The number of reaction and metabolite information of the reconstructed GMN model for the early infection (8h time point) were added to Table 1 in the Results section. The predicted by-product secretion rates based on this metabolic model were also added to Table 2. We also updated Supplementary Figure 8 (PCA of dual transcriptome data) by including the samples of 8th hour of infection.

Submitted filename: Response to reviewers.pdf

Click here for additional data file.

11 May 2022

Dual transcriptome based reconstruction of Salmonella-human integrated metabolic network to screen potential drug targets

PONE-D-21-39848R2

Dear Dr. Cakir,

We’re pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it meets all outstanding technical requirements.

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice for payment will follow shortly after the formal acceptance. To ensure an efficient process, please log into Editorial Manager at http://www.editorialmanager.com/pone/, click the 'Update My Information' link at the top of the page, and double check that your user information is up-to-date. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

Kind regards,

Roman G. Gerlach

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

16 May 2022

PONE-D-21-39848R2

Dual transcriptome based reconstruction of Salmonella-human integrated metabolic network to screen potential drug targets

Dear Dr. Cakir:

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

If we can help with anything else, please email us at plosone@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Roman G. Gerlach

Academic Editor

PLOS ONE