(62c) Flux Balance Analysis at Single-Cell Level

Flux balance analysis (FBA) is a powerful method for simulating cellular metabolism and has been driving genome-scale metabolic network reconstructions. The FBA produces metabolic reaction rates, termed ‘fluxes’, such that these fluxes are stoichiometrically balanced and that they achieve an assumed cellular objective, such as biomass production (Orth, Thiele et al. 2010). FBA has been combined with omics data such as transcriptome (Machado and Herrgård 2014) and proteome (Sánchez, Zhang et al. 2017, Bekiaris and Klamt 2020), to contextualize, constrain, and improve metabolic flux predictions for specific cells. More recently, along with the torrent of single-cell omics data, notably single-cell RNA-sequencing (scRNA-seq) data, there has been a growing interest in applying FBA to simulate metabolic fluxes at single-cell level(Hrovatin, Fischer et al. 2022).

There are challenges in adapting FBA to single cell resolution. First, as it has been shown for bulk data, the expression of metabolic genes is only weakly connected to the flux of reactions that they are associated with based on the gene-protein-reaction (GPR) mapping (Zhang, Li et al. 2010). With the higher noise in scRNA-seq data, this correlation is likely to be lower. Further, scRNA datasets are known to have high sparsity and noise due to technical reasons such as dropouts, as well as biological reasons such as stochastic transcriptional bursting. Therefore, the integration of FBA with scRNA-seq need to account for the complexity of flux control by gene transcription (beyond metabolic genes) and the high level of noise at single cell level.

There have been a few recent developments for flux analysis at single-cell resolution. Yilmaz et al. adapted iMAT, a bulk expression integration method, and developed iMAT++ based on the idea that metabolic fluxes directly correlate with the expression of metabolic genes involved (Yilmaz, Li et al. 2020). Alghamdi et al. developed the algorithm scFEA to generate fluxes of metabolic modules in single cell level using a graph neural network model (Alghamdi, Chang et al. 2021). Here, flux balancing, metabolic gene expression, and cell metabolic activity are implemented in the loss function that is minimized during model training. scFBA is another single-cell FBA method, developed based on linear programming to minimize the sums of a penalty representing the inverse of a gene expression evidence that a reaction is active, while accounting for the stoichiometry within the network (Damiani, Maspero et al. 2019). In these methods, expression of metabolic genes are mapped to the specific reactions based on the GPR association, which ignores the possible roles of non-metabolic genes in (indirectly) regulating metabolic reaction.

In this study, we developed a novel deep learning-based FBA method that leverages the graphical structure of the metabolic network to predict reaction activity from the expression of both metabolic and non-metabolic genes and stoichiometrically-balanced metabolic fluxes. We built on the concept of a reaction flow graph (Beguerisse-Díaz, Bosque et al. 2018) to represent metabolic networks as reaction graphs. Contrary to the typical depiction of metabolic networks where metabolites are nodes and directed edges are reactions occurring between them, in reaction graphs the nodes represent reactions and the edges describe product-substrate relationship. A directed edge between two nodes indicates that the product(s) of the originating reaction (node) is consumed by the destination reaction. The edge weights of the graph capture the probabilistic flow of metabolites between these reactions. The gene expression of a cell is mapped to reaction activity using a single layer perceptron as illustrated in Figure 1A. The discrepancy between the weights in the perceptron and the GPR association is minimized via the loss function. The output of this perceptron is a vector of reaction activity, which serves as the input to the graph neural network layer of the model.

A graph neural network model based on metabolic reaction graph is used to convolve the activity of reactions to predict metabolic fluxes and state (see Figure 1A). By utilizing multiple rounds of graph convolution, the activity of each reaction is passed through the reaction graph using attention-based rule. The number of layers of graph convolution is decided relative to the diameter of the graph (maximum distance between any two nodes in the graph) and through model tuning to allow maximum classification accuracy between the phenotypes. Our rationale is that the flux of a reaction in a cell is determined not only by its noisy gene expression, but also by those of its parent and ancestors. The convolved reaction activity values are used for two tasks. The first is to classify metabolic phenotype that a cell belongs to (e.g. cell type or culture condition). This is to ensure that during the training, the convolved reaction activities are more similar among cells that are of the same or similar phenotype. The second task is to generate stoichiometrically-balanced fluxes. At the single cell level, not all cells may be functioning optimally to satisfy the assumed cellular objective. Thus, a relaxation of flux bounds is done based on the convolved reaction activity. Initially, flux variability analysis (FVA) is performed to identify lower and upper flux bounds that are possible near the optimal flux distribution. During model training, the convolved reaction activities are used to adjust the bounds for reaction fluxes, where reactions with higher (lower) activity are allowed to adopt higher (lower) flux values.

We applied our method to single cell dataset generated by Jariani et al. for Saccharomyces cerevisiae under different growth conditions: glucose vs. maltose media (Jariani, Vermeersch et al. 2020). In this study, yeast cells were initially grown with maltose substrate prior to experiment. At the start of the experiment, it was shifted to a glucose based media for 12 hours. After 12 hours, the substrate was reverted to maltose and grown for another 20 hours. Samples were periodically taken during 5 timepoints and single cell RNA sequencing analysis was made to get single cell transcriptome of the yeast cells. In the study, the authors found that there was a period of time after the substrate was shifted to maltose that growth was arrested. This was a period where the cells were beginning to adapt their metabolism to the new substrate. They termed this period lag-phase. For our study, we utilized the genome-scale metabolic model iMM904 and gene expression data for two time points - glucose at 6 hours(glu-6h) and maltose in lag-phase for 3 hours (Figure 1B) were chosen for determining fluxes. Figure 1B shows a selected set of fluxes which showed notable shifts in fluxes before and after the switch as predicted by our model. We found that hexokinase and ethanol transport and exchange were significantly altered by the switch in media. Indeed, hexokinase (HXK1) gene was reported as being suppressed in the presence of glucose (Jariani, Vermeersch et al. 2020). Further, as were reported by the authors, we do see that biomass flux is nearly the same between the two time points, but the inlag biomass fluxes show a higher variance which indicate that the cells are beginning to adapt to the new substrate and resume growth.

In summary, our multilayered neural network model is able to generate metabolic fluxes that are both stoichiometrically balanced as well as consistent with gene expression at a single cell resolution. Further improvements of this method are currently underway to elucidate novel metabolic signatures for more complex phenotypes and cell subtypes in human cells.

Alghamdi, N., et al. (2021). "A graph neural network model to estimate cell-wise metabolic flux using single-cell RNA-seq data." Genome Research 31(10): 1867-1884.

Beguerisse-Díaz, M., et al. (2018). "Flux-dependent graphs for metabolic networks." npj Systems Biology and Applications 4(1).

Bekiaris, P. S. and S. Klamt (2020). "Automatic construction of metabolic models with enzyme constraints." Bmc Bioinformatics 21(1).

Damiani, C., et al. (2019). "Integration of single-cell RNA-seq data into population models to characterize cancer metabolism." Plos Computational Biology 15(2): e1006733.

Hrovatin, K., et al. (2022). "Toward modeling metabolic state from single-cell transcriptomics." Molecular Metabolism 57: 101396.

Jariani, A., et al. (2020). "A new protocol for single-cell RNA-seq reveals stochastic gene expression during lag phase in budding yeast." eLife 9.

Machado, D. and M. Herrgård (2014). "Systematic Evaluation of Methods for Integration of Transcriptomic Data into Constraint-Based Models of Metabolism." Plos Computational Biology 10(4): e1003580.

Orth, J. D., et al. (2010). "What is flux balance analysis?" Nature Biotechnology 28(3): 245-248.

Sánchez, B. J., et al. (2017). "Improving the phenotype predictions of a yeast genome‐scale metabolic model by incorporating enzymatic constraints." Molecular Systems Biology 13(8): 935.

Yilmaz, L. S., et al. (2020). "Modeling tissue‐relevant Caenorhabditis elegans metabolism at network, pathway, reaction, and metabolite levels." Molecular Systems Biology 16(10).

Zhang, W., et al. (2010). "Integrating multiple ‘omics’ analysis for microbial biology: application and methodologies." Microbiology 156(2): 287-301.