Integrated regulatory and metabolic networks of the tumor microenvironment for therapeutic target prioritization

1 Introduction

Over the past decade, high-throughput single-cell transcriptomics has proven to be instrumental in characterizing within and between cell variation (Tang et al. 2009) and its relation to disease (Ackers-Johnson et al. 2018; Binnewies et al. 2018; Chen et al. 2019; Östman 2012; Paik et al. 2020; Potter 2018), including cancer (Chung et al. 2017; Suvà and Tirosh 2019; van Galen et al. 2019; Zhu et al. 2017). The tumor microenvironment is a diverse, complex collection of many interacting cell types (Whiteside 2008). Such diversity necessitates highly individualized therapies. The push toward personalized medicine has been recognized as a critical step toward effective cancer treatments (Chin et al. 2011; Tannock and Hickman 2016). Despite the progress and availability of high-throughput technologies, translating genomic discovery to clinical endpoints has not kept pace and remains a longstanding bottleneck in the field (Chin et al. 2011). This is partly due to the associated time and expense of experiments needed for drug development. Drug discovery is time-consuming and costly; recent estimates to bring a single cancer drug to market are $648 million and revenues are more than three times that (Prasad and Mailankody 2017).

Computational systems biology can be used to accelerate the discovery process through the integration of computational and experimental research (Kitano 2002a,b). Computational experiments, known as in silico experiments, are used to make predictions, generate hypotheses, and prioritize the most promising experiments and treatments for future investigations. Network analysis is a cornerstone of systems biology that leverages mathematical models of biological networks to enhance our understanding of function, mechanisms, and disease. Structures of these network models are often obtained either from experimental studies, learned from the data, or some combination thereof (Liu et al. 2020). Networks typically target specific domains of biological variables such as protein-protein interaction (Kovács et al. 2019), metabolism (Angione 2019), or gene regulation (Ferreira et al. 2018). This work focuses on two classes of networks that model different biological domains: regulatory signaling and metabolic networks, and addresses their integration. The mathematical frameworks used for modeling these domains are fundamentally different, which has made their integration a longstanding challenge.

Ordinary differential equations are used to model cellular metabolism. These models depend on mathematical representations of biochemical reactions governed by mass balance equations with a well-understood stoichiometry (Nelson et al. 2008). Dynamic systems require kinetic rate parameters that are often unknown and difficult or impossible to measure. Furthermore, cost and challenges of omics data collection at multiple time points have limited its use in the post-genome era. A popular alternative is to model the system at steady-state, where there are no changes in metabolite concentrations, and the model reduces to a linear system. Steady-state models can be parameterized using observational data from a single time-point. Constraint-based modeling approaches are used for flux analysis in steady-state systems (Gudmundsson and Thiele 2010; Kauffman et al. 2003; Orth et al. 2010; Palsson 2015). The models rely on an objective function, such as biomass production, together with thermodynamic constraints to define region of feasible fluxes for optimization. Constraint-based models have become widespread and accessible due to the large-scale published computational models of cellular metabolism used for metabolic modeling. The models consist of metabolic reactions in a particular organism catalyzed by one or several enzymes encoded by genes.

Constraint-based models have been used to characterize differences in flux distributions between malignant and normal cells, thereby providing insights into cancer-specific metabolic hallmarks and enzyme targets for anti-cancer treatment (Agren et al. 2014; Folger et al. 2011; Frezza et al. 2011; Gatto et al. 2015; Yizhak et al. 2014; Zhao et al. 2013). Notably, the adaptation of these models to different cancers is justified because the differences between the mRNA level observed under different steady-state conditions have been shown to explain most of the variation in concentrations of the associated proteins (Li and Biggin 2015; Liu et al. 2016b). Although the encoding is often not directly modeled, there is a direct gene-protein-reaction association assumption that enables the use of high-throughput genomic data within the frame of metabolic networks (Bordel 2018; Gu et al. 2019).

Regulatory pathways and processes have been curated in databases (Ashburner et al. 2000; Kanehisa and Goto 2000). However, unlike metabolic stoichiometry, structures of regulatory pathways are often less understood. Inference of gene regulatory networks using graphical modeling paradigms has been a very active area of research (Bansal et al. 2007; Benfey and Mitchell-Olds 2008; Friedman 2004; Husmeier 2003b; Karlebach and Shamir 2008; Lee and Tzou 2009; Markowetz and Spang 2007; Menéndez et al. 2010; Rockman 2008; Vijesh et al. 2013). The bulk of research in this area is dedicated to inferring the network structure from transcriptomics data to understand direct or indirect associations between genes. New methods like Functional Inference of Gene Regulation (FigR) use single-cell ATAC-seq and RNA-seq data to infer gene-regulatory networks, examining chromatin dynamics and gene expression for a unique view of regulatory networks (Kartha et al. 2022). This work utilizes a directed class of Probabilistic Graphical Models (PGMs) known as Bayesian Networks (BNs) (Koller and Friedman 2009; Lauritzen 1996) for models of gene regulation. A valuable asset within the BN paradigm is the use of probabilistic reasoning, which enables inferences about probabilistic changes within the network under different conditions and perturbations (Moharil et al. 2016; Yu et al. 2020).

There have been limited attempts to link transcriptomics data into metabolic models. Context-specific metabolic models incorporate information from transcriptomics data into flux estimation. The underlying premise is that not all reactions are uniformly essential or active across different cell types or conditions, and that transcriptomics carries this information as a measure of abundance. Boolean models have been used to embed transcriptional regulation in the constraints (Covert and Palsson 2002; Covert et al. 2008; Shlomi et al. 2007). More complex models directly use transcriptomics data to constrain the bounds of estimated fluxes (Chandrasekaran and Price 2010; Marmiesse et al. 2015), generating multi-valued logical rules for genetic regulation with steady-state FBA models.

These constraint-based models of metabolism that utilize transcriptomics fail to account for the complex structure of the regulatory networks that influence the enzymes facilitators of biochemical reactions. Ignoring the hierarchical structure of the regulatory network could indirectly invoke an independence assumption between enzymes in the model. To this end, there has been significant progress in the mathematical modeling of regulatory networks (Brazhnik et al. 2002; Dondelinger et al. 2013; Hill et al. 2012; Huynh-Thu et al. 2010; Schäfer and Strimmer 2005; Smolen et al. 2000; Zhang and Horvath 2005), which has been largely untapped in metabolic models. BNs have been widely used to model regulatory networks, primarily focusing on inferring structures (Blair et al. 2012; Hageman et al. 2011; Husmeier 2003a,b; Wang et al. 2013; Werhli and Husmeier 2007; Werhli et al. 2006). Recently, Yu and Blair (2019) developed a novel framework that utilizes BN models to perform exact probabilistic reasoning using bulk transcriptomics (Koller and Friedman 2009). This framework recasts the network as a computational model for perturbation and prediction. This framework was used to predict brain metabolism changes after perturbations to the HIF-1A pathway in Alzheimer’s patients.

The integration of single-cell transcriptomics within computational models of metabolism has been limited. Damiani et al. (2019) proposed a method for single-cell flux balance analysis (scFBA) to infer single-cell fluxomes. The approach utilizes single-cell transcriptomics, environmental constraints relating to medium-cell flux exchanges, and assumptions on cell interactions. Earlier work by Damiani et al. (2017) developed population FBA (popFBA), which utilized scRNA to characterize the heterogeneity in the tumor environment for classic FBA models. To the authors’ knowledge, there has not been progress in integrating scRNA regulatory networks with models of metabolism.

This study is motivated by the observation that drug interventions targeting transcription factors (TFs) impact the metabolism differently across cells in the tumor microenvironment (DeBerardinis 2020; Reina-Campos et al. 2017). This can be partially explained by differential cell signaling among cell types (Arneth 2020; Balkwill et al. 2012; Kenny et al. 2007). Accounting for these differences will provide mathematical models that are more physiological. These models can be used to advance our understanding of the differential effects of therapeutic perturbations on metabolisms in different cell types in the tumor microenvironment.

This work develops a modeling paradigm for prediction in integrated regulatory and metabolic networks. The bridge between the two paradigms is probabilistic reasoning in BNs, which generates predictions to constrain the metabolic model. These cell-specific integrated models provide a unique opportunity to characterize the metabolism of different cell types in response to upstream genetic perturbations, providing insights on therapeutic target prioritization, formulation of combination therapies, and future drug discovery. This model integration framework is generalizable to other applications, such as different cell types, organisms, and diseases.

Applications focus on networks parameterized using single-cell RNA sequencing data from different cell types in the glioblastoma tumor microenvironment (Neftel et al. 2019) and a metabolic model of the brain. The model centers on the HIF-1A pathway, which is central to cancer metabolism (Chiavarina et al. 2010; Koukourakis et al. 2002; Rawłuszko-Wieczorek et al. 2014) and includes a complex network of transcription factors, many of which are drug targets. In silico experiments are performed to examine the effects of 177 curated TF-targeting drugs published in the drug repurposing hub library (Corsello et al. 2017) on metabolism of cells in the tumor microenvironment. Results show that drug combinations inhibiting STAT3 (e.g., Niclosamide) with other transcription factors (e.g., AR inhibitor Enzalutamide) are predicted to strongly suppress anaerobic metabolism in malignant cells, without major interference to other cell types’ metabolism, suggesting a potential combination therapy for anti-cancer treatment.

2 Materials and methods

The integration of regulatory and metabolic networks is challenging due to the fundamental differences between the two modeling frameworks. Section 2.1 describes the modeling of the regulatory network using probabilistic PGMs. Probabilistic reasoning in BNs is described, which will enable predictions of downstream enzymatic changes after perturbation. Section 2.2 describes constraint-based models of metabolism used to predict the fluxomes. Sections 2.3 and 2.4 describe model integration and an application to cancer, respectively. Lastly, a sensitivity analysis is presented in Section 2.5, to assess the sensitivity of predictions to network construction and the integration of model systems.

2.4 Model integration and pharmaceutical interventions evaluation

Cell-specific BN models were integrated with a metabolic model to evaluate the system-wide effects of pharmaceutical interventions targeting TFs. The model is used to infer the impact of TF perturbations on the regulatory network and downstream metabolism. The integration of these two models is achieved by embedding predicted fold changes of enzyme levels from the signaling pathways into the constraints of the metabolic model for flux variability analysis. The integration follows closely to Yu and Blair (2019) but differs in terms of the models themselves and the intervention strategy.

Figure 1 shows a schematic of this integrated model and process. As described in Section 2.3.2, cell-specific BNs are inferred based on prior knowledge about the network structure (Garcia-Alonso et al. 2019; Holland et al. 2020) and participation in the metabolic model (Figure 1A) (Gavai et al. 2015). These BNs are perturbed using belief propagation to alter the TF levels. The TF perturbations and combinations of perturbations were guided by drug mechanisms of action (MOA) data of 6798 drugs from the drug repurposing hub (DRH) library (Corsello et al. 2017). The MOA data recorded the specific biological mechanism (e.g., reaction or interaction) that produces the pharmacological effect. A total of 177 inhibitors were mapped such that the BN models cover all their known targets. Supplementary (S-Table 1) details each drug’s target TFs and MOA. Pharmaceutical interventions of all 177 curated drugs were evaluated using the cell-specific BNs.

Figure 1:

A schematic describing the integrated modeling approach. (A) Cell-specific models are developed using Bayesian Network inference using prior constraints of regulatory structure from a database and single-cell RNA data. The system-wide effects of perturbations to these models are explored using probabilistic reasoning for BNs. The enzyme levels estimated from these models are embedded into the constraints for flux variability analysis. (B) Lsei-FBA is performed using a metabolic model of brain metabolism. The outputs are predicted cell-specific reaction flux rates that can be compared to estimated baseline fluxes, which provides insights into the relative changes that occur in the cell-specific systems under TF perturbations.

2.4.1 Cell-specific enzyme fold change

For each cell-specific network, belief propagation was performed to predict the expression level of the enzymes and TFs in the regulatory network before and after the TF perturbation (Figure 1B) (Lauritzen 1996). We employed the BP procedure in the BayesNetBP package (Yu et al. 2020), which implements the belief propagation algorithms described in (Cowell 2005). The output of interest is cell-specific enzyme fold changes for the constraint of the flux variability analysis, which were calculated based on the following three-step procedure.

1. Initial enzyme expression level: The imputed expression data were utilized to infer both the structure and the parameters of cell-specific BNs. Cluster trees for belief propagation were then inferred for each network. The initial expression levels of enzymes and TFs were estimated as the marginal means (Yu et al. 2020).

2. Mutilated network after TF perturbation: Pharmaceutical interventions of each drug were treated as the perturbations (knockdown) of the target TFs. To simulate the effect of drug interventions, we assume the activity of all drug target TFs will be reduced by 2 standard deviations from their basal (average) level. Moreover, in order to model the causal relationship after drug interventions, a mutilated graph was obtained by removing all inward connections to the TF nodes that are targeted by the intervention (White and Vignes 2019). Cell-specific drug-intervened cluster trees were then generated using the mutilated graphs for propagation.

3. Enzyme expression level after TF perturbation: The perturbed TF activities after drug intervention were treated as evidence and absorbed into the cluster trees from the mutilated BNs. Belief propagation was performed for the estimation of perturbed enzyme levels. Cell-specific enzyme fold changes were then calculated using the mean expression level after perturbation (numerator) relative to the initial enzyme level (denominator).

2.4.2 Cell-specific relative flux change

Cell-specific relative changes for the metabolic fluxes were calculated in order to assess the effects of pharmaceutical interventions targeting TFs on metabolism. Calculations follow the three-step procedure described below.

1. Basal flux estimation in oligodendrocytes: Baseline flux rates in oligodendrocytes were estimated through the maximization of net ATP production, which reflects the high requirements for energy production in the brain, which is critical for bioenergetics, function, and neurodegeneration (Du et al. 2008). This objective function was also applied in other studies of brain metabolic modeling (Gavai et al. 2015; Yu and Blair 2019). Following (Gavai et al. 2015; Yu and Blair 2019), the flux in oligodendrocytes before the intervention using FBA with total ATP production maximized and with the pentose phosphate pathway (PPP) and GABA shunt fluxes constrained to measured values relative to glucose uptake. Basal flux in oligodendrocytes will then be used as the input for the Lsei-FBA model to estimate the flux distributions in oligodendrocytes and other cell types before and after interventions.

2. Enzyme fold change and initial cell-specific flux estimation: Cell-specific enzyme levels for all four cell types before and after TF perturbations were obtained following the same BP procedure described in Section 2.4.1, resulting in eight total estimates. Each of these enzyme estimates was contrasted with baseline oligodendrocytes to obtain fold change estimates. For each cell type, initial estimates for the fluxes were calculated as the product of basal flux in the oligodendrocytes for that reaction and the enzyme fold change in the associated cell-specific enzyme expression. In cases where multiple enzymes mapped to a single reaction, the average fold change across these genes was used to constrain the corresponding flux.

3. Estimation of cell-specific relative flux change: Employing the Lsei-FBA model in the R LIM package (Haskell and Hanson 1981; Van Oevelen et al. 2010), the final prediction of cell-specific metabolic fluxes was obtained by minimizing the sum of squared errors between estimated and predicted flux (basal flux × fold-change of enzymes). The constraints were (1) the steady-state assumption for intercellular metabolites that are not exchanged between brain and blood, (2) the restriction to forward flux through irreversible reactions (inequality constraints). The cell-specific relative flux change (C) is calculated as:

   C = v 1 − v 0 v 0 ,

   where v ₀ and v ₁ are the fluxes at baseline and after perturbation, respectively. Relative flux change of 0 means no change, while a negative value indicates that the drug down-regulates the reaction.

3 Results

BNs were learned for each of the four cell types using the scRNA-seq data and prior graph constraints generated from databases. Unconnected upstream TF nodes were dropped within each cell-specific network. The BNs of four cell types showed distinct learned structures that support cell-specific cell-signaling. The macrophage network (Figure 2A) had 94 nodes and 174 edges; the malignant cell network (Figure 2B) had 103 nodes and 207 edges; the oligodendrocyte network (Figure 2C) had 81 nodes and 129 edges; the T-cell network (Figure 2D) had 76 nodes and 104 edges. The learned BNs of malignant cells and macrophages were larger in terms of nodes and edge density, suggesting a more intricate signaling network when compared to oligodendrocytes and T-cells.

Figure 2:

The inferred BNs by cell-type: (A) macrophage, (B) malignant, (C) oligodendrocyte and (D) T-cell. Enzymes in the metabolic model are indicated in red.

Measures between graphs provide additional insights into the similarities between cell-type structures (Figure 3). Notably, since the BNs were learned under the same structural prior, larger networks with more edge density are likely to be more similar in their structure. BNs of malignant cells and macrophages had the most shared edges (N = 98), while sparser networks of oligodendrocytes and T-cells had the least (N = 28) shared edges (Figure 3A). The Jaccard index computes shared edges adjusting for the total edge number of the union of two graphs. Therefore, it provides a normalized measure of similarity between networks (higher is more similar). Jaccard indices showed similar results as shared edges that BNs of malignant cells and macrophages were the most similar (Jaccard index = 0.35), while oligodendrocytes and T-cells were the most dissimilar (Jaccard index = 0.14). In addition, the greatest difference was found between BNs of malignant cells and T-cells with a Jaccard index of 0.15 (Figure 3B). Graph edit distance (GED) which measures the number of changes (i.e., adding, removing, substituting nodes and edges) to transfer one graph to the other one, was also utilized to capture the similarity between networks (lower is more similar). Malignant cells and T-cells were the most different, with the largest GED of 229, while oligodendrocytes and T-cells were more similar, with the least GED of 177 (Figure 3C). Finally, closeness centrality distance (lower is more similar) captures the difference between two networks with respect to a given node centrality, thereby accounting for the local structures around different nodes in the networks. In gene regulatory networks, a gene with high centrality is often viewed as a hub with important biological importance. Overall, we found that malignant cells are the most different from the others, with the highest closeness centrality distances from other cell types (Figure 3D), indicating the mutated transcriptional dysregulation in malignant cells.

Figure 3:

Heatmaps of similarity and dissimilarity measurements between four inferred cell-specific BNs: (A) number of shared edges (higher is more similar), (B) Jaccard index (higher is more similar), (C) graph edit distance (lower is more similar) and (D) closeness centrality distance (lower is more similar). Similarity and dissimilarity measurements are indicated in red and blue, respectively.

Systematic perturbations were made to the integrated model in order to explore the system-wide effects of interventions that target TFs on the regulatory network and downstream metabolism. In this study, the primary focus is the enzymes and flux rates from three pathways: glycolysis, anaerobic glycolysis, and pyruvate dehydrogenase (PDH). Since each network is uniquely parameterized by scRNA-seq data, we observed unique belief propagation results, despite structural similarities. Patterns of the predicted changes of downstream enzymes (Figure 4) and fluxes (Figure 5) were visualized after TF perturbations. For visualization purposes, the set of perturbations was selected based on the magnitude and diversity of fold changes. The full set of 247 perturbations from target TFs, and combinations of target TFs from 177 drugs can be found in Supplementary (S-Table 1, S-Figure 1 and S-Figure 2).

Figure 4:

Clustering of the predicted fold changes in enzymes when transcription factors and combinations of transcription factors in the network are targeted (rows). There were five distinct clusters observed (A–E). The fold changes are systematically displayed for the enzymes for the malignant, T-cell, macrophage and oligodendrocyte (columns). This ordering is preserved across the different pathways: glycolysis, anaerobic glycolysis and the flux from the glycolytic pathways to the TCA cycle (pyruvate dehydrogenase).

Figure 5:

Clustering of the predicted relative changes in fluxes when transcription factors and combinations of transcription factors in the network are targeted (rows). There were five distinct clusters observed (A–E). The flux changes were aggregated and are systematically displayed for the malignant, T-cell, macrophage and oligodendrocyte (columns). This ordering is preserved across the different pathways: glycolysis, anaerobic glycolysis and the flux from the glycolytic pathways to the TCA cycle (pyruvate dehydrogenase).

Large variations were observed that indicate a strong differential effect across both TFs and cell types (Figure 4). Clustering the differential effects of the predicted enzyme expression after TF target inhibition revealed five distinct groups (Figure 4A–E). Inhibition of TFs promoting the response to hypoxia is expected to divert the flux from lactate production to oxidative phosphorylation, which is consistent with our findings. However, the ways that different cells achieving were observed to be markedly different. Patterns of inhibition were observed for combinations of NFKB1 with other TFs (Figure 4A). Enzymes from all three reaction groups: glycolysis, anaerobic glycolysis and PDH increased in macrophages, T-cells, and were most pronounced in oligodendrocytes. On the other hand, malignant cells exhibited modest reductions in desired anaerobic glycolysis and PDH with combination therapies of SP-100030 (JUN and NFKB1 inhibitor) and drugs that target on AR (e.g. Enzalutamide, DHEA sulfate, Epitiostanol, Bicalutamide, etc.).

Drugs in cluster B (Figure 4B) mainly targeted on JUN together with other TFs. Notably, oligodendrocytes were not sensitive to JUN inhibition, while glycolysis enzymes such as PFKP and PGK1 were upregulated only in malignant cells. Cluster C (Figure 4C) showed that inhibition of RELA with other TFs leads to significant upregulation in LDHA and PDHB in only oligodendrocytes. Moreover, targeting ESR1 with other TFs (TP53, KLF4 or RARA) caused enzymes decrease only in malignant cells, while not much effect on the other cells. Cluster D (Figure 4D) includes the drug combinations that target on HIF1A (e.g. Acriflavine, Bakuchiol, etc.) with other TFs. Under these drug interventions, enzymes in glycolysis and anaerobic glycolysis showed major decreases in malignant and macrophage cells, while T-cells and oligodendrocytes were not as sensitive. Inhibition of STAT3 (e.g. Niclosamide, Meisoindigo, Napabucasin, etc.) with other TFs were clustered as group E (Figure 4E), and showed similar but weaker effect patterns as inhibitions of HIF1A (Figure 4D).

Relative flux rate changes under combined drug interventions also showed distinct patterns (Figure 5). For simplicity, the predicted flux changes within glycolysis, anaerobic glycolysis and PDH were averaged within cell-type and the effects of perturbations. Clustering analysis yielded five distinct groups (Figure 5A–E), which is consistent with predicted enzyme fold changes (Figure 4). Cluster A corresponds to drugs targeting RELA or NFKB1 with other TFs (Figure 5A) that caused an increase in the glycolysis and PDH reactions, and a decrease in anaerobic glycolysis only in oligodendrocytes cells, while other cells were not found to be as sensitive. Cluster B showed that inhibition of targets combinations that led to anaerobic glycolysis decrease in all four cell types (Figure 4B). Patterns of inhibition of ESR1 with select TFs (FOXM1, RARA, KLF4, TP53) are shown in Cluster C (Figure 5C), which only down-regulated anaerobic glycolysis in malignant cells, but not in other cell types. STAT3 inhibition in combination with other TFs (Figure 5D) formed cluster D, which shows decreases in anaerobic glycolysis in malignant cells, and much more slightly in the other three cell types. Under inhibition of HIF1A with other TFs in cluster E (Figure 5E), reactions across all pathways were also down regulated in both malignant and macrophage cells, and anaerobic glycolysis was decreased across all cell types. The sensitivity analysis affirmed the consistency of these findings (S-Figure 3). Particularly, STAT3 inhibition alongside other TFs results in decreased anaerobic glycolysis in malignant cells, and slightly less so in the other three cell types. This robust pattern validates the effectiveness of our proposed integrated model.

Taken together, the most promising targets arise from drug combinations inhibiting STAT3 with other TFs due to the strong effects on the malignant cells, and minimal effects on other cell types. To further investigate the targets, we selected a promising drug combination from STAT3 cluster (Figures 4E and 5D) that corresponds to the combined drug intervention of Enzalutamide (targets AR) and Niclosamide (targets STAT3). Figure 6 visualizes the estimated relative flux rate change of the four cell types under this combined intervention. Relative change of 0 means no change, and negative relative change indicates down-regulation. Results showed that the drug intervention has almost no effects on the PDH reactions, and it only slightly restrained glycolysis in macrophages and malignant cells with relative flux changes of −0.15 and −0.14, respectively. However, this combined treatment strongly suppressed anaerobic reaction in malignant cells with a relative change of −0.92, while the down-regulation effects for other cell types were much milder, with flux fold change ranging from −0.14 to −0.34. Taken together, these results suggest that Enzalutamide and Niclosamide drug combinations may be a potential treatment that is specific in suppressing anaerobic metabolism in malignant cells, without major interference to other cell types.

Figure 6:

The inferred relative flux changes for each of the four cell-types after inhibition of Enzalutamide (targets AR) and Niclosamide (targets STAT3).

Our finding that the STAT3 inhibitors have specific effects on glioblastoma cells is consistent with the public drug sensitivity data and literature. The primary PRISM repurposing dataset (Corsello et al. 2020) obtained from the cancer dependency map in the DepMap portal (Ghandi et al. 2019) shows that most glioblastoma cell lines are sensitive to the STAT3 inhibitors (S-Figure 4). Recent literature also found that napabucasin, a STAT3 inhibitor, can suppress the proliferation and invasion of glioblastoma cells (Han et al. 2019). On the other hand, hyperactivation of STAT3 in tumor-infiltrating immune cells causes immunosuppression by inhibiting both innate and adaptive immune responses (Zou et al. 2020).

4 Discussion

This work develops a novel approach for integrating cell-specific gene networks into constraint-based models of metabolism. To the authors’ knowledge, this represents the first approach to integrating genetic cell-specific networks based on single-cell transcriptomics data into constraint-based models of metabolism. This integration represents a more physiological system and allows us to characterize the distinct metabolic profiles within the tumor microenvironment under different perturbations. In the context of cancer, these perturbations can be viewed as potential single drug targets or combination therapies. Although this work focused on cancer, the general framework can be adopted in other applications, such as different cell types, model organisms, and diseases.

A strength of this work is the flexible BN representation of signaling networks to estimate the system-wide impacts of perturbations. There are several advantages to using BNs in this context. For structural learning, prior information about regulatory associations can be used to constrain the learning of the graph structure. Probabilistic reasoning utilized exact inferences that are only possible because of the directed acyclic nature of the graphs (Lauritzen 1996). Although the directed nature of the graphs is a strong assumption, BN structures have been shown to mimic biological mechanisms and capture feedback (Blair et al. 2012). Moreover, exact reasoning is not possible for undirected networks.

On the other hand, a unique challenge for BN learning and inference is the zero inflation of the scRNA-seq data (Cao et al. 2021; Luecken and Theis 2019), which violates the Gaussian distribution assumptions for the parameterization of CGBNs. The zero-inflation issue poses a universal challenge to basic analytics performed with scRNA-seq and has been shown to have biological relevance (Choi et al. 2020). In this work, we proposed a hybrid imputation method leveraging enzyme expression level and TF activities that overcome the issues related to inference with zero-inflated data. We applied MAGIC algorithm (Van Dijk et al. 2018) to impute and smooth the enzyme expression levels and VIPER technique (Alvarez et al. 2016) to infer the TF activities based on the TF transcript abundance. There are two benefits in using inferred TF activities. First, the TF protein activity is a more direct predictor for enzyme levels than the TF transcript abundance. Second, the distribution of TF activities can be better modeled by Gaussian distributions. A future direction will be directly utilizing the zero-inflated data, which would require alternative parametric and modifications for the inference.

The integrated computational model was used to explore a perturbation space of TFs. Pyruvate dehydrogenase (PDH) is a gatekeeper enzyme between glycolysis and the TCA cycle. The modulation of PDH expression has been actively investigated as a target for controlling cancer growth (Sradhanjali and Reddy 2018). In order to prioritize our findings, we sought to identify target perturbations that inhibit PDH in malignant cells and minimize the dysregulation of the metabolism of other cells in the tumor microenvironment. Some of these promising findings link to drugs in current use.

The androgen receptor (AR) is a steroid hormone activated TF that is targeted by Enzalutamide. In normal tissues, androgen effects are mediated by AR through hormone-receptor binding in both males and females. Androgen-independent AR activation is a common finding in castration-resistant prostate cancers. Enzalutamide has been used primarily for prostate cancer (Beer et al. 2014; Tran et al. 2009). AR has also been considered as a potential target for glioblastoma (Yu et al. 2015; Zalcman et al. 2018), the most common form of brain cancer. The incidence of glioblastoma is 1.6 times higher in men than women (Carrano et al. 2021), and the upregulation of AR has been implicated in promoting tumorigenesis (Yu et al. 2015). Despite these associations, there has only been a recent study that confirms the therapeutic role of targeting AR in glioblastoma in vitro and in vivo using Enzalutamide in mouse models (Zhao et al. 2021).

STAT3 is central in regulating the anti-tumor immune response and plays essential roles in inhibiting the expression of crucial immune activation regulators and promoting the production of immunosuppressive factors (Zou et al. 2020). Niclosamide was the only launched STAT3 targeting drug in our drug list and is known to have broad clinical applications for the treatment of various cancers (Chen et al. 2018). A study showed that Niclosamide synergistically reversed Enzalutamide resistance in prostate cancer cells. Combined treatment with Niclosamide plus Enzalutamide resulted in inhibition of colony formation and growth arrest by inducing cell apoptosis (Liu et al. 2015). Multiple studies investigated the effect of Niclosamide in combination therapy. For example, Niclosamide with other drugs like terpyridine or cisplatin achieved significant treatment effects in breast cancer (Karakas et al. 2015; Liu et al. 2016a). Niclosamide combined with erlotinib was shown to be effective in treating head and neck cancer (Li et al. 2013a) and lung cancer (Li et al. 2013b).

This approach also has limitations. In the modeling, the scope of the metabolic model is limited to a simple representation of central carbon metabolism. This work represents a tractable proof of concept. Genome-scale metabolic models have been developed and made available (Schellenberger et al. 2010). A future direction of this work will be to expand upon the integrated modeling framework to ensure scalability. Flux balance analysis also requires an objective function, which in our applications was to maximize ATP production for basal flux estimation in oligodendrocytes, reflecting the high energy demands of brain metabolism. Setting an objective function is more straightforward in single-cell systems where the objective is growth. In human cells, a single objective is unrealistic. Importantly, when simulating the effect of drug perturbations, the current approach does not consider potential off-target effects. Drug actions can often influence other targets beyond the intended ones, which can cause unanticipated changes in cellular functions. Model validation is challenging in integrated systems. These potential limitation should be kept in mind when interpreting our findings and could be an additional aspect to consider in future model refinements.

In conclusion, we have developed a flexible framework for integrating cell-specific BNs of signaling pathways into metabolic models using scRNA-seq data. This generalizable approach leverages BN inference to make predictions about the system under perturbations, and flux estimation to explore the effects of perturbations on the metabolic system. Metabolic models and gene regulatory networks have both been widely studied, but in parallel. Biologically, their integration is well understood, but mathematically connecting these models has been a longstanding challenge. Network integration is a promising approach for bridging these domains into a computational model that can be used to make predictions and gain novel insights about disease.