Integration of Metabolic Modeling with Gene Co-expression Reveals Transcriptionally Programmed Reactions Explaining Robustness in Mycobacterium tuberculosis

Abstract

Robustness of metabolic networks is accomplished by gene regulation, modularity, re-routing of metabolites and plasticity. Here, we probed robustness against perturbations of biochemical reactions of M. tuberculosis in the form of predicting compensatory trends. In order to investigate the transcriptional programming of genes associated with correlated fluxes, we integrated with gene co-expression network. Knock down of the reactions NADH2r and ATPS responsible for producing the hub metabolites, and Central carbon metabolism had the highest proportion of their associated genes under transcriptional co-expression with genes of their flux correlated reactions. Reciprocal gene expression correlations were observed among compensatory routes, fresh activation of alternative routes and in the multi-copy genes of Cysteine synthase and of Phosphate transporter. Knock down of 46 reactions caused the activation of Isocitrate lyase or Malate synthase or both reactions, which are central to the persistent state of M. tuberculosis. A total of 30 new freshly activated routes including Cytochrome c oxidase, Lactate dehydrogenase, and Glycine cleavage system were predicted, which could be responsible for switching into dormant or persistent state. Thus, our integrated approach of exploring transcriptional programming of flux correlated reactions has the potential to unravel features of system architecture conferring robustness.

Figure 1. Identifying correlated reaction sets in Genome Scale Metabolic Reconstruction.

(A,B) The genome scale metabolic network reconstruction updated iNJ661m (iNJ661mu) built using genome annotation, and bibliographic data. Using stoichiometric coefficient as thermodynamic constraint and other constraints as reaction bounds, these models could be used to simulate growth rate or biomass production (objective function in LP problem) in silico using Flux Balance Analysis approach assuming steady state. (C) Represents the wild type flux distribution of reactions. (D) Matrix of flux profiles of reactions shows flux profile of reaction R1 and flux profiles of affected reaction R2. All fluxes are shown V1 and V2 in WT and V1a is flux of R1 at 2% knock down and V2a is new flux of R2 at 2% flux reduction of WT flux of R1 and so on. Example is given as IGPS (Indole-3-glycerol phosphate synthase) knock down and affected reaction PIt (Inorganic phosphate transporter).

Figure 2. Illustrative example of positively and negatively correlated reactions in constraint based chemical reaction model.

X – axis: Flux knockdown of flux of R1, and Y-axis: Flux values. In the upper panel the flux profile of knocked down reaction is shown. Lower panel shows two types of reactions. In positively correlated reactions (R2, R3 and R4) with reduced fluxes on knock down of R1, these constitute the fragility in the system. In negatively correlated reactions (R5, R6 and R7) with increased fluxes on knock down of reaction (R1), these constitute the robustness in the system. A real example of Cysteine synthase (CYSS) knock down and positively and negatively affected reactions is shown with the schematic representation.

Figure 3. Numbers of reactions with correlated fluxes corresponding to knocked down of reactions in different metabolic pathways.

X-axis: the pathways to which the knocked down reaction belongs and Y-axis: Numbers of reactions of affected pathways shown as bar histograms (Blue color: Positively correlated reactions and Red color: Negatively correlated reactions).

Figure 4

2-dimensional colourmap display of the scale of effects scored in terms of the numbers of reactions affected (reaction pairs) with respect to the total number of possible pairs in affected pathways on knockdown of reactions of a given pathway (in percentages) (A) positively correlated reactions, (B) negatively correlated reactions. The colour code corresponding to the numbers of affected reactions is shown in the corresponding scale bar. The knockdown pathways are in the rows and the affected pathways are in the columns.

Figure 5. Gene expression patterns of subunits of enzymes (ATP Synthase, NADH dehydrogenase and Mycobactin synthase) encoded by multiple genes in 520 microarray samples.

(X-axis = microarray data samples, and Y-axis = Z-score expression values) (A) The subunits of ATP synthase are encoded by 8 genes. All these 8 genes are co-expressed and have positive correlation in expression values. (B) The NADH dehydrogenase is encoded by 14 genes and all 14 genes are co-expressed and have positive correlation in expression values. (C) The Mycobactin synthase is encoded by 7 genes and all 7 genes are co-expressed and have positive correlation in expression values.

Figure 6. Work flow of integration of gene co-expression with reaction correlation.

Figure 7. Proportion of reaction pairs under co-expression in metabolic pathways.

X-axis: pathways and Y-axis: ratio of reaction pairs with associated genes co-expressed to the total number of reaction pairs. Blue colored bar: positively correlated reactions, Red colored bar: negatively correlated reactions.

Figure 8. Gene expression correlation of genes associated with reactions of glycolysis and citric acid cycle pathways.

The correlation values of one gene versus the rest are plotted in color coded format. The color code range bar corresponding to the correlation values is shown. (A) fumarase (B) malate dehydrogenase (C) glyceraldehyde-3-phosphate dehydrogenase (D) Phosphofructokinase (E) phosphoglycerate kinase (F) triose phosphate isomerase. Rows: genes of knocked down reaction, Columns: genes associated with negatively correlated flux of reactions with associated genes co-expressed. The scale bars are uniformly set −1 to +1 in colour scale ranging from red to green in all cases. However, in cases where the extreme values are absent, the colours are only assigned to the next minimum and maximum values available. Therefore the interpretations is that red signifies normalised low score values whereas green signifies normalised high score values.

Figure 9. Reciprocal gene expression between the duplicated genes coding for inorganic phosphate transporter (PIt).

(A) Patterns of gene expression between the two PIt genes Rv0545c and Rv2281 in 520 samples. X-axis: Microarray sample, Y-axis: Z-score of the gene expression values. (B) Heatmap of PCC between PIt associated genes and genes associated with negatively correlated flux reactions during Pit knockdown. Rows: PIt genes, Columns: genes associated with negatively correlated reactions.

Figure 10

Reciprocal gene expression between the duplicated genes cysteine synthase (CYSS) (A) Gene expression patterns of three copies of cysteine synthase gene. The gene Rv0848 is negatively correlated in expression with Rv2334. The genes Rv2334 and Rv1336 had positive correlation in expression. X-axis: Microarray sample, Y-axis: Z-score of the gene expression values. (B) Heatmap of PCC between CYSS associated genes and genes associated with negatively correlated flux reactions during CYSS knockdown. Rows: CYSS genes, Columns: genes associated with negatively correlated reactions.