article: Metabolic reprogramming is a general feature of cancer cells. Regrettably, the comprehensive quantification of metabolites in biological specimens does not promptly translate into knowledge on the utilization of metabolic pathways. By estimating fluxes across metabolic pathways, computational models hold the promise to bridge this gap between data and biological functionality. These models currently portray the average behavior of cell populations however, masking the inherent heterogeneity that is part and parcel of tumorigenesis as much as drug resistance. To remove this limitation, we propose single-cell Flux Balance Analysis (scFBA) as a computational framework to translate single-cell transcriptomes into single-cell fluxomes. We show that the integration of single-cell RNA-seq profiles of cells derived from lung adenocarcinoma and breast cancer patients into a multi-scale stoichiometric model of a cancer cell population: significantly 1) reduces the space of feasible single-cell fluxomes; 2) allows to identify clusters of cells with different growth rates within the population; 3) points out the possible metabolic interactions among cells via exchange of metabolites. The scFBA suite of MATLAB functions is available at https: //github. com/BIMIB-DISCo/scFBA, as well as the case study datasets.
  Cancer is a heterogeneous, multi-factorial and essentially genetic disease, in which various types of mutations alter the functioning and interactions of genes, causing cancer cells to proliferate in an uncontrolled manner. Despite the plethora of cancer-related mutations, a reduced number of recognizable phenotypic hallmarks [1,2] have been identified. Metabolic rewiring, in particular, is a general feature of cancer cells, which reprogram their metabolism to feed their unrestrained proliferation, as it requires high amounts of energy and building-blocks [3]. The design principles underlying the causative role of metabolism in promoting growth as a function of the nutritional constraints are starting to be investigated [4,5] and the idea of targeting the distinctive features of cancer metabolism has received considerable attention [6]. Unfortunately, a single metabolic program cannot be used to globally define an altered tumour metabolism [7], as cancer cells, even within the same tumour, may cope with the above metabolic requirements by engaging different metabolic pathways [8]. Such variability produces different dependencies on exogenous nutrients, and reflects into heterogeneous responses to metabolic inhibitors [9]. Furthermore, in solid tumours, cancer cells are embedded within the tumour microenvironment (TME), a complex network of fibroblasts, myofibroblasts, myoepithelial cells, vascular endothelial cells, cells of the immune system and extracellular matrix. TME also includes chemical gradients of oxygen and nutrients: the complex interaction of all these elements plays a major role in tumour metabolic heterogeneity [10]. The metabolic interplay that occurs among cancer cells and TME—supported by experimental evidence on how malignant cells may extract high-energy metabolites (e. g. , lactate and fatty acids) from adjacent cells [11,12]—contributes to treatment resistance [13]. Therefore, effective therapeutic strategy should incorporate knowledge of intra-tumour metabolic heterogeneity and cooperation phenomena within cancer cell populations. Knowledge about the utilization of metabolic pathways requires quantification of metabolic fluxes (i. e. , the rate at which a substance is transformed into another through a given reaction or pathway). As quantification of the full complement of cell metabolites (metabolomics) alone does not provide information on internal fluxes [14], they can be measured only indirectly, mainly via several isotope-labeled metabolomics experiments coupled with metabolic flux analysis. A more functional option is Flux Balance Analysis (FBA), which uses linear algebra algorithms to solve a mass balance problem, given constraints on the flux of some relevant reactions, with particular regard to extracellular fluxes (rate of intake and secretion of metabolites) [15]. An advantage of FBA is that, as opposed to intracellular fluxes, extracellular fluxes can also be approximated from measurements of the concentration of metabolites in the spent cell culture medium at different time points. FBA has proved able to correctly predict the growth yield of microorganisms when extracellular fluxes are constrained [16–18]. Several approaches have been proposed to set further constraints on internal fluxes, by exploiting other -omics data, such as transcriptomics or proteomics data. Protein abundances should be used cautiously as proxy for flux rates, as the availability of substrates ultimately determine the reaction rate. Transcript levels should be used even more prudently than protein abundances, because many factors beyond transcript concentration contribute to determine the expression level of a protein, such as translation and degradation rates, spatial locations and post-transcriptional regulation. However, the differences between the mRNA level observed in different conditions (at the steady state) explain most of the variation in concentration of the associated proteins [19,20]. Therefore, coupling the information of transcripts with that on extracellular fluxes, within the FBA steady-state modeling framework, provides a reasonable solution to the problem of predicting intracellular fluxes. Nevertheless, extracellular fluxes can be hardly measured at the single-cell level. They cannot be approximated from the spent cell culture medium, because it portrays the extracellular fluxes of the bulk. They might be measured with isotope-labeling and metabolic flux analysis techniques, but this would require single-cell metabolomics. Unfortunately, single-cell metabolomics is still at its embryonic stage, mainly because of limitations in working with minute amounts of material [21–23]. This is a major problem, because flux distributions are currently estimated from metabolic measurements retrieved from bulk samples, which often contain intermixed and heterogeneous cell subpopulations, thus overlooking possible cooperation and compensation phenomena. In the simplest example, let us imagine that two populations exist: one that secrete lactate and another that consumes all the lactate produced by the former. The predicted flux distribution will represent the sum of the two populations, hence the flux through lactate production/consumption will be regarded as inexistent, portraying a behavior far from the real one. This limitation might partially be overcome by analyzing a subpopulation of cells supposedly having a more homogeneous metabolism, for example, by using fluorescence-activated cell sorting techniques to isolate it according to specific physical properties. However, when dealing with populations of cells derived from tumours, it is difficult to assess the relative composition of intermingled cancer cells and of stromal elements within the tumor architecture, and fluorigenic markers may not robustly correlate with metabolic phenotype. Hence, stratification, in terms of metabolic function and of consequent susceptibility to metabolic drugs, of heterogeneous populations taken, for example, from biopsies, xenografts or organoids, requires characterization of cellular metabolism at the single-cell level. At this purpose, we here introduce the first computational framework scFBA (single-cell Flux Balance Analysis) to predict single-cell fluxomes and possible metabolic interactions among them, starting from (bulk) extracellular fluxes and single-cell transcriptomes. As compared to single-cell metabolomics, the field of single-cell transcriptomics has indeed progressed to a deeper resolution [5,24]. An established classification [25–29] of current methods to integrate (bulk) transcriptomic data into constraint-based models is based on the purpose of the method, and distinguishes approaches that aim at extracting the active metabolic network (typically from a genome-scale one, e. g. , [30–33]) from methods that aim at estimating the extent of the flux of each reaction in a network (e. g. , [34–39]). However, these methods cannot be directly extended to single-cell modeling, because of the lack of information of extracellular fluxes at the single-cell level discussed above. Even if it was possible to measure single-cell extracellular fluxes, there is no protocol to measure transcriptome and fluxes of the very same single-cell (the cell is destroyed after either analyzing the trascriptome or metabolome). For this reason, we propose to use a multi-scale model, in order to solve a unique mass balance problem to identify the possible combination of single-cell steady states that concurrently satisfies constraints on single-cell transcriptomes and extracellular bulk fluxes. In order to determine the flux distribution of each network, each cell is not considered in isolation: it is allowed to interact with other cells in the population, via release/uptake of metabolites into/from the TME. Datasets obtained from cancer biopsies or patient-derived xenografts are the ideal candidate to capture the heterogeneous composition of tumors. As currently no information on scRNA-seq and extracellular fluxes on the very same sample is publicly available, as a proof of principle, we here applied scFBA to lung adenocarcinoma (LUAD) patient derive xenograft scRNA-seq, collected by [40], while testing different sets of constraints for the extracellular fluxes. To prove the robustness and applicability of our method, we also applied scFBA on further independent breast cancer datasets collected by [41].
  Although some attempts to study the cooperation between different metabolic populations have been put forward [42–44], mostly focused on microbial communities, these methods require indeed a priori knowledge about the specific metabolic requirements and objectives of the intermixed populations. Unfortunately, even though metabolic growth may approximate the metabolic function of some cell populations, we cannot assume that each cell within an in vivo cancer population proliferates at the same rate, nor that it proliferates at all. A major example is given by the different proliferation rates of stem and differentiated cells [45]. For this reason, differently from other approaches [44], we do not impose that the population dynamics is at steady-state (and hence that cells all grow at the same rate), although we do continue to assume that the metabolism of each cell is. Conversely, scFBA aims at portraying a snapshot of the single-cell (steady-state) metabolic phenotypes within an (evolving) cell population at a given moment, and at identifying metabolic subpopulations, without a priori knowledge, by relying on unsupervised integration of scRNA-seq data. We have previously shown how Flux Balance Analysis of a population of metabolic networks (popFBA) [46] can in line of principle capture the interactions between heterogeneous individual metabolic flux distributions that are consistent with an expected average metabolic behavior at the population level [46]. However, the average flux distribution of a heterogeneous population can result from a large number of combinations of individual ones, hence the solution to the problem of identifying the actual population composition is undetermined. To reduce this number as much as possible, we here propose to exploit the information on single-cell transcriptomes, derived from single-cell RNA sequencing (scRNA-seq), to add constraints on the single-cell fluxes. An identical copy of the stoichiometry of the metabolic network of the pathways involved in cancer metabolism is first considered for each single-cell in the bulk. To set constraints on the fluxes of the individual networks, represented by the single-cell compartments of the multi-scale model, we took inspiration from bulk data integration approaches that aim to improve metabolic flux predictions, without creating context-specific models from generic ones [34–39]. At the implementation level, we use continuous data, rather than discrete levels, to overcome the problem of selecting arbitrary cutoff thresholds. At this purpose, some methods (e. g. [30,32]) use expression data to identify a flux distribution that maximizes the flux through highly expressed reactions, while minimizing the flux through poorly expressed reactions. To limit the problem of returning a flux distribution (or a content-specific model) that does not allow to achieve sustained metabolic growth, we use instead the “pipe capacity” philosophy embraced by other methods, such as the E-Flux method [36,37], of setting the flux boundaries as a function of the expression state. These methods tend to use relative rather than absolute expression values. For instance, the original formulation of E-flux [36] sets relative boundaries in relation to the most expressed reactions. In order to avoid comparing enzymes with different gene-protein translation rates, which may also largely differ in their kinetic parameters (e. g. binding affinity) and in the number of associated isoforms/subunits, we prefer to normalize boundaries in relation to the condition/cell/tissue in which a given reaction is mostly expressed, as done in a more recent version of the E-flux method [37] and in other continuous methods [34,35,47]. Specifically, we distribute the total (bulk) possible flux of each reaction proportionally to the activity score of that reaction in each cell. To compute a score for reactions that involve many genes, similarly to other approaches [36–38,48], we assume that enzyme isoforms contribute additively to the overall activity of a given reaction, whereas enzyme subunits limit its activity, by requiring all the components to be present for the reaction to occur. Alternative but similar approaches, such as [30], consider the maximum value (instead of the sum) of isoform values. A scheme of the scFBA approach is depicted in Fig 1. Here we briefly recall the popFBA approach. For a more comprehensive description, the reader is referred to [46]. Starting from a template metabolic network map A, corresponding to a generic single-cell and defined as A = (X A, R A, E A) —where X A = { X 1, …, X M } is the set of metabolites in network A, R A = { R 1, …, R N } the set of biochemical reactions taking place among them and E A = { E 1, …, E N ext } is a set of Next unbalanced reactions (exchange reactions), enabling a predefined set of metabolites (including the pseudo-metabolite representing biomass) Y = { Y 1, …, Y N ext } ⊂ X A to be inserted in or removed from the system—the popFBA procedure first builds a population model composed of Ncells replicates Ac of network A, each one corresponding to a single-cell c, c = 1, …, Ncells, and which can cooperate by exchanging nutrients in the tumour microenvironment. For each single-cell c, A c = (X c, R c, C c) is its metabolic network, where: Accordingly, a new set of metabolites pertaining to the TME compartment Y ′ = { Y i ′ } with i = 1, …, Next must also be defined. Because original exchange reactions have been replaced by cooperation reactions, a new set of Nblood exchange reactions B = { B 1, …, B N blood } is defined, which allows a subset of metabolites K = { K 1, …, K N blood } ⊂ Y ′ to be exchanged with the external environment, e. g. , the blood supply: The population model P is then defined by: (i) the union set of the metabolites X P = ∪ c X c ∪ Y ′; (ii) the internal reactions R P = ∪ c R c; (iii) the cooperation reactions C P = ∪ c C c; (iv) the population exchange reactions B. A stochiometric matrix SP is then built for all reactions in R P, C P and B and for all metabolites in X P and Y ′. The final size of matrix SP is (Ncells ⋅ M + Nblood) × (Ncells ⋅ (N + Next) + Nblood). Once the population model is obtained, the total biomass of the Ncells single-cells is maximised by means of linear programming, as in standard FBA [15]. The solution of popFBA represents the flux distribution v → = (v 1, …, v N cells · (N + N ext) + N blood) = (r 1 1, …, r N 1, …, r 1 N cells, …, r N N cells, c 1 1, …, c N ext 1, …, c 1 N cells, …, c N ext N cells, b 1, …, b N blood) that maximises the biomass exchange flux bbiomass, with vi representing any flux i of the population model, and for each single-cell c, r i c representing the i-th internal flux, c i c representing the i-th cooperation flux and bi an exchange flux with blood. The optimization problem is postulated as follows: maximise b b i o m a s s subject to S v → = 0 →, v → L ≤ v → ≤ v → U (3) v → L and v → U are vectors specifying the lower and upper bound respectively for each flux vi of v →. A negative lower bound indicates that flux is allowed in the backward reaction. To solve the above problem we exploited the Gurobi solver within the COBRA Toolbox [49]. scFBA takes as input a template metabolic network map A, as in popFBA, plus a scRNA-seq dataset in the form of a Ngenes × Ncells matrix T, where Ngenes is the number of genes and Ncells is the number of single-cells under study. Each element Tg, c, g = 1, …, Ngenes, c = 1, …, Ncells corresponds to the normalized read count of gene g in cell c such as, for instance, the TPM (Transcripts Per Kilobase Million). The risk of the presence of false negatives in RNA-seq, and in particular scRNA-seq, is an established problem. Although a totally safe solution does not exist, scFBA allows to employ the information on bulk expression profile, when available, to manage the risk, by envisioning the following scenarios. We define a Reaction Activity Score (RAS), for each single-cell c = 1, …, Ncells, and each reaction j ∈ R, based on Gene-Protein-Reaction association rules (GPRs). GPRs are logical formulas that describe how gene products concur to catalyze a given reaction. Such formulas include AND and OR logical operators. AND rules are employed when distinct genes encode different subunits of the same enzyme, i. e. , all the subunits are necessary for the reaction to occur. OR rules describe the scenario in which distinct genes encode isoforms of the same enzyme, i. e. , either isoform is sufficient to catalyze the reaction. In order to compute the RAS we distinguish: In case of composite reactions, we respect the standard precedence of the two operators. The first step of the scFBA approach is the creation of a multi-scale population model, composed of Ncells, according to the popFBA described above, but starting from the template metabolic network A*. As described in the subsection related to input and data pre-processing, A* is a subnetwork of the generic model A which integrates the transcriptional information that holds for all cells in the bulk. Once the population model is obtained, the scFBA approach imposes two kinds of constraints: In order to project the information of the activity score of a given reaction j, in a given cell c, R j c, onto its flux “pipe capacity”: In this work, we mainly use the following 3 LUAD datasets obtained from the NCBI Gene Expression Omnibus (GEO) data repository under accession number GSE69405. We repeated all the analyses on the following independent breast cancer datasets (GEO access number: GSE75688), including scRNA-seq data of single-cell suspensions of cancer tissues obtained on the day of the surgery of untreated breast cancer patients [41]: Each of the 5 datasets includes the gene expression level of more than 20. 000 genes in the form of Transcript Per Kilobase Milion (TPM). We filtered out a few cells with less than 5000 genes detected. For each dataset, we retained only the metabolic genes included in HMRcore model (418 genes). The dataset transcripts are identified by Ensembl ID, which we automatically converted into HUGO Gene Nomenclature Committee (HGNC) ID. The datasets also contain the expression profile of the bulk samples, which we used to pre-process data as described above. From the computational perspective, the scFBA approach is suitable for simulation of genome-wide metabolic networks, such as [51,52]. However, in view of previous analyses [4], in order to have more control on the analyses and make the interpretation of results more straightforward, at this stage, we preferred to focus on a more handful and carefully reconstructed core metabolic network. We used, as template network A, the metabolic core model HMRcore introduced in [53] and used in [46,48]. As exchange of fatty acids between cells in tumours has been recently reported [12,54], we included the possibility to exchange palmitate via the TME and, accordingly, mitochondrial palmitate degradation and gluconeogenesis. Given the importance of reactive oxygen species (ROS) metabolism observed in [4], we also inserted ROS production and removal pathways. As the original version of the model does not include information on GPRs, such rules have been extracted from Recon 2. 2 [51] and included in the HMRcore model. We decided to disregard the GPRs associated to the complexes I to IV of the electron transport chain in scFBA computations, because it unrealistically requires up to 81 genes (AND rule). However the flux through complexes I to IV should be modulated by the constraints on complex V (ATP synthase). The final version of the HMRcore model includes 315 reactions (of which 263 are associated with a GPR) and 418 metabolic genes. The SBML of the model is provided in https: //github. com/BIMIB-DISCo/scFBA.
  We first applied scFBA to the 5 datasets described in the Methods section, assuming maximization of total (population) biomass synthesis rate as objective function. All five population models displayed a non negligible maximal growth rate, something that cannot be taken for granted when integrating transcriptomics into FBA models [25,27,29]. In S1 Fig and accompanying S1 Text, we also show that, if ϵ takes value 0, the scFBA problem is still feasible, but we obtain very small values for the fluxes. In order to highlight that the scFBA approach efficiently reduces the space of optimal solutions, we compared the variability of the biomass production flux of each of the Ncells single-cells simulated within the population model, for each of the 5 datasets under study. We report in Fig 2 the results for the two datasets relative to the primary tumors, and in S2 Fig for the other datasets. When no information on cells’ transcriptome is employed (as in standard popFBA settings [46]), the type 2 constraints of the metabolic network are identical for all cells. This implies that each cell is capable of contributing alone to 100% (100) of the objective function value (i. e. , the biomass of the total population). As depicted in Fig 2 (left plots) and S2 Fig (left plots), the biomass flux value of each cell, within the set of optimal solutions, spans indeed from 0 to 100% of the total biomass (purple rectangles). On the contrary, after scRNA-seq data integration, as performed via scFBA, the biomass flux of each cell can only take a specific (optimal) value, corresponding to a certain fraction of the total biomass (green rectangles, which results in a single line, because the maximum and minimum optimal flux values coincide.) To show how this volume reduction in the space of alternative optima may actually affect predictions, we performed a single gene deletion analysis with and without scRNA-seq data integration (scFBA and popFBA, respectively). When a single gene is deleted, the reactions for which the expression of such gene is essential (i. e. , reactions exclusively associated to the gene, or reactions associated to that gene and other genes with an AND operator) are removed from the network (i. e. , from the set R c, ∀ c). After removal, the population model is newly optimized for total biomass production, and the growth ratio (grRatio) of the new biomass over the previous one is computed. Fig 2 and S2 Fig (right bar plots) report the grRatio observed for those genes deletions that displayed a different effect before and after data integration. Notice that when the grRatio equals 0, the corresponding bar is not displayed at all. To verify that the differences between scFBA and popFBA are not a mere consequence of the removal of reactions (in scFBA) that are inactive in all cells of the bulk from the template metabolic network A, we include in the plots the prediction of the isolated template metabolic network A*. We refer to this third simulation setting as bulkFBA. However, bulkFBA includes information on on-off reactions only. It is not possible to modulate the flux capacity of reactions as a function of gene expression, because it not possible to compute relative expression values. Remarkably, some genes that are redundant (grRatio = 1) in popFBA settings (i. e. , with no scRNA-seq data integration) may even become essential in scFBA settings (i. e. , with scRNA-seq data integration) (grRatio = 0). This is the case, in lung adenocarcinoma, of the following genes: HGNC: 10980, which encodes enzymes responsible for glutathione/phosphate, fumarate/phosphate or α-ketoglutarate/malate antiports; HGNC: 8806, which encodes for a subunit of pyruvate dehydrogenase; HGNC: 8896, which encodes for an isoform of phosphoglycerate kinase and HGNC: 4141, which encodes for an isoform of glyceraldehyde 3-phosphate dehydrogenase. In breast cancer, only gene HGNC: 8806 falls into this category. Conversely, some genes that display a significant effect (grRatio ≈ 0. 5) in popFBA become instead redundant in scFBA. This is the case, in lung adenocarcinoma, of the genes that encode for ATP synthase (HGNC: 823,830,833,837,838,840,845-848,14247,8016), indicating that the integration of scRNA-seq data forces a (suboptimal) flux distribution for cancer cells which, consistently with the well-known Warburg effect, does not rely largely on ATP synthase for ATP production, thus resulting in a milder effect when the reaction is depleted. Worth of note, although these genes are not completely redundant (grRatio < 1) in breast cancer, the deletion of the respiratory chain has indeed a mild effect in both tumors, as well as in the other cancer datasets reported in S2 Fig. It is apparent, from Fig 2 and S2 Fig (right bar plots), that bulkFBA provides intermediate results between scFBA and popFBA. Some genes that are redundant in PopFBA, are lethal in both scFBA and bulkFBA. This is the case for example in LCPT45 of gene HGNC: 10980 (Mitochondrial dicarboxylate carrier). This result is expected, given that its isoform has been deleted according to bulk data in both simulations. Conversely, some gene deletions that have a significant effect according to bulk data have no effect when also single-cell data are considered, in particular the genes encoding for components of the respiratory chain. Also, the effect of the deletion of pyruvate kinase (HGNC: 9021) is smaller in scFBA than in bulkFBA. On the other hand, some deletions may show some effect only when scRNA-seq are considered. This is particularly true for genes that are involved in cooperation mechanisms among cells, as for instance gene HGNC: 29, in both datasets, whose product promotes the secretion of palmitate, which can be taken up by other cells. As previously mentioned, single-cell fluxes are expected to be less noisy than transcript signals, which are typically analyzed by means of multi-variate statistical analysis [5] and, therefore, the former might be used to better identify cell clusters that might represent distinct metabolic subpopulations. To confirm this hypothesis, we performed a cluster analysis on the expression values (scRNA-seq) of the metabolic genes and compared the results with those of a cluster analysis performed on the fluxes predicted by scFBA. To this end, we performed both hierarchical and k-means cluster analysis. In order to avoid reactions with typical high flux-value, or genes with high expression, to induce a bias on clustering results, we first remapped the flux (transcript) values of each reaction (gene) j in the interval [0,1]: value 0 is assigned to the cell showing the lowest value for a given flux (transcript), 1 to the one showing the highest value. Fig 3 and S3 Fig report the results of the hierarchical clustering analysis (distance metric: euclidean), for transcripts (left column) and fluxes (middle column), respectively for the two primary tumors and for the other 3 datasets under study. From the dendrograms and heat-maps, one can see that cells cluster in a few well-separated groups, when the extracted features (the fluxes) are considered, whereas they cluster more in “singletons”, when the original features (the transcripts) are used. For example, when observing the fluxes computed for the LUAD dataset LCPT45 (panel A in Fig 3), it is apparent that two major (and a minor) groups of cells can be identified, corresponding respectively to the blue and red-coloured leaves in the dendrogram. We evaluated which are the most different fluxes among the two major groups, by using the Z score test of statistical significance (S1 Table). Remarkably, the two groups significantly differ in their growth rates (Z-scores: 3. 2). 82 reactions significantly differ between the two groups with a 99% confidence level. This set mostly include pathways directly or directly linked with biomass synthesis, such as biosynthesis of fatty, amino and nucleic acids. To quantitatively compare the clustering of transcripts and fluxes, we first performed a k-means clustering with different number of clusters k, by considering n = 100 bootstrap iterations (with random centroid assignments) and by selecting the clustering resulting in the maximum inter-cluster distance. We then assessed the clustering goodness, by means of traditional “elbow” and “silhouette” evaluation methods. We refer to S2 Text for details about these approaches. The elbow method in the right column of Fig 3, indicates that a elbow is observed at k = 3 for the fluxes relative to the primary tumour datasets, hence the optimal number of clusters is 3, which corresponds indeed to the k identified by the hierarchical clustering analysis for these two datasets. In S3 Fig (panel D), we evaluated the silhouette for the dataset LCPT45 transcripts (left) and fluxes (right) for k = 3, i. e. , the value identified from the elbow analysis, which also corresponds to the highest average silhouette value, when varying k in {2, ⋯, 6}. Consistently with the more ready recognition of major clusters in the flux diagrams noted above, the drop in the sum of squared errors (SSE) is much stronger and the average silhouette value is considerably higher in the flux case than in the transcripts case (where the average coefficient is close to 0) indicating that the calculation of fluxes leads to a better clustering as compared to the evaluation of transcripts. All in all, the results of the cluster analyses indicate that fluxes can be better clustered than their transcript counterparts. Remarkably, it can also be noticed that the LUAD primary tumour xenograft (see the sharper “elbow” and the clustergram in Fig 3A) fluxes better partition into clusters than the fluxes of the cell line (S3 Fig, H358) and of the secondary tumour xenograft (S3 Fig, LCMBT15). This result is in line with the data reported in [40], indicating that a binary separation of cells is more evident in LCPT45. We detected a similar difference between the clustering results of primary and secondary tumour of the independent breast cancer datasets (BC04 in Fig 3 and BC03LN in S3 Fig). Indeed, the former population is more heterogeneous in the binary sense than the latter [41]. The main rationale behind solving a unique mass balance problem for many cells together, given constraints on the extracellular fluxes of the bulk, rather than many separate mass balance problems, is that the nutrient consumption and secretion rates (extracellular fluxes) can be readily determined or approximated from measurements of the concentration of metabolites in the cell culture media at different time points for the bulk only. Another major side benefit of this approach is that it allows to identify the possible interactions among cells within a population, as pointed out in [46]. We verified that, after data integration, some cells secret metabolites that are up-taken by other cells. The heat map in Fig 4 shows the (normalized) flux values of cooperation reactions for the LCPT45 dataset: a positive value means that the cell is secreting the metabolite in the tumour microenvironment, whereas a negative flux that the cell is uptaking it from the tumour microenvironment. It can be observed that a complex network of interactions is established among cells. In particular, a consistent group of cells consumes the lactate and palmitate that are secreted by other groups. The scatter plots in Fig 4 show the dispersion of the fluxes of uptake/secretion from/into the TME for lactate and palmitate and how they couple with different growth rates, portraying a relationship far more complex than that depicted with popFBA (no scRNA-seq integration [46] and no exchange of palmitate allowed). More in detail, the dominant subpopulation (larger cluster), which includes the majority of the cells, displays secretion of NH3 and glutammate and uptake of lactate and palmitate. Within the dominant subpopulation, only mild differences in flux rates can be detected, whereas within the other (smaller) subpopulations heterogeneous patterns of nutrient production/consumption are observed, which result from metabolic interactions. It is also worth remarking that the predictions obtained with scFBA would not be possible with context-specific extraction methods from bulk data. In this respect, we compared the results of scFBA with two widely-used methods, GIMME [62] and iMAT [63], by focusing on the extracellular fluxes of the dominant single cell population identified with scFBA (see S2 Table for details). As expected, GIMME and iMAT cannot predict the consumption of certain metabolites, unless they are provided as exogenous nutrients. Indeed, both methods would not predict the consumption of lactate and palmitate. Even the prediction on produced metabolites may differ from that of the dominant population in scFBA: although both methods predict NH3 to be produced, consistently with scFBA predictions, glutamate would be produced according to iMAT only. Metabolic interactions between cancer-associated fibroblasts and cancer cells, mediated by palmitate [12,54] and lactate [59–61] have been recently reported. At the same time, it has been shown that metabolic heterogeneity can arise in genetically homogeneous cells as simple as the budding yeast Saccharomyces cerevisiae [64]. scFBA has the potential to highlight possible metabolic heterogeneity also within a genetically homogeneous population of cancer cells. The validation of the predicted interactions requires however non-trivial ad hoc experiments. As current techniques do not allow for easy determination of metabolites at the single-cell level, the heterogeneous population of alive cells should be first sorted to separate it into the sub-populations identified by scFBA. However, to sort cells based on fluorescent labeling (Fluorescence activated cell sorting), further analyses are necessary to possibly identify markers differentially expressed by the sub-populations. Less direct approaches might be taken, for example, by measuring the growth rate of wild-type cell populations and mutant populations in which the cooperation has been prevented (e. g. by blocking secretion or uptake of involved metabolites) and comparing with the model predictions. At this purpose, it should be assessed whether the metabolic interactions identified by scFBA are actually advantageous for tumor growth or are just related to the entropy of the system, i. e. , to the fact that a configuration in which interactions among heterogeneous phenotypes take place is more likely than a configuration of identical and independent phenotypes. We address the issue in the following paragraphs. Both in popFBA and in scFBA the cells were able to collaborate metabolically. As the integration of scRNA-seq data greatly reduced the space of feasible FBA solutions, those data encode information on how nutrient utilization should be distributed amongst the individual cells. Some cells that can no longer carry out a certain part of a pathway let their neighbors do this. However, it should still matter which nutrients are available to all cells. For a deeper characterization of given cancer populations, exo-metabolomic measurements to constrain the population boundary conditions would thus be needed. An exhaustive sensitivity analysis of scFBA results to boundary conditions is out of the scope of this work. However, it is interesting to compare the conditions in which the two major metabolites involved in cooperation (i. e. , lactate and palmitate) are externally supplied to the population or must be produced endogenously. Notably, we observed that uptake of exogenous palmitate does not affect the biomass production rate, indicating that no growth advantage is conferred by free availability of lipids. This result is in line with experimental evidence that cancer cells rely on de novo synthesis of palmitate-derived lipids [3]. However, in the baseline setting (no external palmitate supplied), we observed a group of cells that uptake the palmitate secreted by others (Fig 4). We verified that, once internalized in those cells, palmitate is not processed by the beta-oxidation pathway, but directly contributes to the biomass synthesis, supporting the evidence reported in [68], that an exogenous source of fatty acids can substitute for de novo synthesis in promoting cell proliferation and attenuate the cancer-specific toxic effect of lipogenesis inhibitors. It has also recently been pointed out that a limited access to environmental lipids may render the cancer cells more sensitive to the inhibitors of lipogenesis [68]. In line with these findings, it can be observed, with regard to the LCPT45 population (Fig 5A), that a set of genes stops being essential when exogenous palmitate is supplied. As expected, this set mainly includes genes directly involved in the synthesis of palmitate, namely: citrate synthase, fatty acid synthase and pyruvate dehydrogenase. The expansion of the plot in Fig 5A (left) shows that the latter (pyruvate dehydrogenase, ID: HGNC: 8808) is essential for each cell within the population. It should indeed be noted that when the synthesis of palmitate is prevented in all cells, and exogenous palmitate is not supplied, also cells that used to rely on the palmitate synthesized by other cells must be affected. As opposed to palmitate, the metabolite lactate is not strictly required for growth. However, it can be observed in Fig 5B that the deletions of genes encoding for glucose-6-phosphate isomerase (HGNC: 4458) and for triosephosphate isomerase 1 (HGNC: 12009) —two important steps for the utilization of glucose through glycolysis are not essential when lactate uptake is allowed, suggesting that lactate may be able to replace glucose as carbon source. Interestingly, when lactate uptake is prevented, the plot expansion in Fig 5B (left) shows that the gene HGNC: 4458 is essential in many but not all cells. Also the set of metabolites allowed to be released, e. g. , in the blood may affect the effect of gene deletions. For instance, if glutamate secretion is prevented, the deletion of the gene that encodes palmitate secretion becomes essential, as shown in Fig 5C. Remarkably, it has been reported that secretion of lipids facilitates tumour progression [69], whereas inhibitors of glutamate release have been proposed as new targets for breast cancer-induced bone-pain [55]. scFBA may enable to shed light on how the disposal of carbons through these two metabolites relates with the utilization pattern of exogenous nutrients.
  We have here introduced scFBA to solve the problem of reconstructing the potential single-cell fluxome, starting from single-cell transcriptomes, by taking into account environmental constraints, as well as cell-cell interactions. Importantly, scFBA is able to point out the metabolic interactions that are established within a cell population. scFBA integrates sc-transcriptomics data with (bulk) extracellular fluxes of the same cancer cell population, by means of a computational approach inspired to complex systems science [70]. A limitation of our approach is that it uses mRNA levels as a proxy of the maximal velocity of reactions (Vmax), thus neglecting the many factors that contribute to determine the expression level of a protein [71], as well as the role played by binding affinities of proteins in determining the Vmax. However scFBA does not predict the single-cell fluxome as a linear function of the assumed Vmax: the constraints on mass-balance and on the availability of substrates, as determined by the rates of consumption/secretion of nutrients by the entire population and by the requirement for the tumour mass to grow, are simultaneously taken into account. scFBA might also be implemented by using sc-proteomics rather than transcriptomics, should the former become available. Although we do not explicitly model spatial organization, the constraints on single-cell transcriptome should implicitly preserve the information on the usage/secretion of most nutrients of each cell in its original position. The method seems also to neglect communication between cells through growth factors. In reality, it mostly does not if the growth factors act by changing transcription. But if they act by phosphorylating enzymes, then this is not taken into account. In its current form scFBA can however already be implemented to investigate why metabolic drugs are often ineffective and to provide indications for more effective treatment. As a proof of principle, we have successfully applied the methodology to LUAD and breast cancer datasets. We have shown that the integration of scRNA-seq greatly reduces the space of feasible solutions that sustain metabolic growth of the overall population, which is prerequisite of tumour growth. This reduction allowed us to restrict the set of candidate drug targets, by eliminating targets that may seem obvious for the bulk, but do not work for heterogeneous populations, and, on the other hand, by revealing targets whose relevance can be appreciated only if cooperation among heterogeneous cells is accounted for. We have also illustrated that scFBA is valuable to extract features (i. e. , the sc-flux values) from scRNA-seq data, in order to identify metabolic clusters of cells, which may be used to investigate other fingerprints of the cancer metabolic deregulation. Although popFBA assumes that the cells achieve optimal growth, this assumption is mitigated in scFBA, by taking into account the transcriptional constraints. Moreover, we have shown how alternative objective functions, such as ATP maximization, may be investigated. Sub-optimality may also be taken into account by using sampling methods [72]. In this study, we have used scRNA-seq obtained with protocols based on C1 Single-Cell Auto Prep System, which have the advantage of allowing to remove dead cells before sequencing but may suffer from including low numbers of cells per sample (34—55 cells), as compared to modern emerging technologies which allow to obtain the single-cell transcriptome of thousands of cells at the cost of a bulk experiment [73], with an improved number of genes/transcripts per cell. In the future, our approach may be readily generalized to this kind of data. As illustrated in S4 Fig, the time of a scFBA computation increases linearly with the number of simulated cells and with the size of the template metabolic network. Alternatively, bootstrap-like methods can be adopted to reduce the number of cells considered in a single computation, while parallelizing the simulation of many smaller systems. A major challenge, when dealing with scRNA-seq, is the presence of false-negatives. As a first approximation, we have used a bulk RNA-seq expression filter, where single-cell expression values of transcripts never detected in scRNA-seq but detected in bulk RNA-seq are replaced by the bulk values. For more reliable genes, we preserve information on cell heterogeneity, but we mitigate the risk of false-negatives setting the bound of the associated reactions to a small value ϵ rather than completely removing it. The choice of the value of ϵ is partially arbitrary. However, we verified that main results of our work are robust with respect to this choice (S1 Text and S1 Fig). The methodology might be refined, by combining it with more sophisticated data pre-preprocessing techniques, which may also take into account the specific quality parameters of the dataset. Finally we have shown how constraints on the nutrient consumption and secretion rates (extracellular fluxes) of the specific sequenced population may affect scFBA predictions. As opposed to intracellular fluxes, (bulk) extracellular fluxes might be readily estimated, e. g. , by approximation from metabolite concentrations in spent medium (exo-metabolome), when culturing patient-derived cells. Hence, measurements of both extracellular fluxes and single-cell transcriptional information of the same heterogeneous cancer populations are needed to make scFBA predictions fully reliable. These datasets can be realistically obtained by culturing population of cancer cells and then processing the cells with scRNA-seq technologies and analyzing their spent medium with biochemical analyzers. Experiments under a controlled setting, e. g. , a co-culture of metabolically characterized cancer cell lines, may first be performed to validate and tune the capability of scFBA to identify and measure the prevalence of different metabolic subpopulations of cells. The application of scFBA to analyze datasets more representative of tumor heterogeneity, obtained, for instance, by culturing cells from human biopsies, xenografts or organoids, will then pave the way to cancer personalized medicine.
summary: Cytotoxicity of chemotherapeutic agents and resistance to targeted treatments are the main reasons why cancer is still one of the top causes of death. As tumor cells are intrinsically resistant to therapies that target signaling pathways, targeting the metabolic hallmarks of cancer holds promise for more incisive treatments. Regrettably, the heterogeneity of cancer metabolism hinders the identification of effective treatments. To fully uncover the metabolic heterogeneity within tumors, characterization of metabolic programs (metabolic flux distributions) at the single-cell level is required. To fill the gap between current technologies for genomics and future technologies for fluxomics, both at the single-cell and the genome-wide scale, we propose to integrate cancer data from: 1) single-cell transcriptomics and 2) bulk metabolomics, into a multi-scale stoichiometric model, to deliver for the first time metabolic fluxomes at the single-cell level. To this end, we introduce a new paradigm for flux balance analysis and data integration in cancer metabolism to: 1) characterize metabolic heterogeneity, not only at the inter-, but also at the intra-tumor level 2) identify the metabolic interactions between cancer populations, whose role in resistance to metabolic treatments has been recently recognized 3) predict the collective response to drug targeting of metabolism.
keywords: cell physiology
  medicine and health sciences
  metabolic networks
  cell metabolism
  physiological processes
  oxygen metabolism
  metabolites
  network analysis
  genome analysis
  pharmacology
  drug metabolism
  computer and information sciences
  metabolic pathways
  pharmacokinetics
  biochemistry
  cell biology
  physiology
  transcriptome analysis
  genetics
  biology and life sciences
  secretion
  genomics
  metabolism
  computational biology
title: Integration of single-cell RNA-seq data into population models to characterize cancer metabolism