Optimal control nodes in disease-perturbed networks as targets for combination therapy

Abstract

Most combination therapies are developed based on targets of existing drugs, which only represent a small portion of the human proteome. We introduce a network controllability-based method, OptiCon, for de novo identification of synergistic regulators as candidates for combination therapy. These regulators jointly exert maximal control over deregulated genes but minimal control over unperturbed genes in a disease. Using data from three cancer types, we show that 68% of predicted regulators are either known drug targets or have a critical role in cancer development. Predicted regulators are depleted for known proteins associated with side effects. Predicted synergy is supported by disease-specific and clinically relevant synthetic lethal interactions and experimental validation. A significant portion of genes regulated by synergistic regulators participate in dense interactions between co-regulated subnetworks and contribute to therapy resistance. OptiCon represents a general framework for systemic and de novo identification of synergistic regulators underlying a cellular state transition.

Fig. 1

Network controllability theory and overview of the OptiCon algorithm. a Identification of driver nodes in a directed graph based on structural controllability theory. By converting a directed graph G into a bipartite graph, a structural control configuration (SCC) of the network can be identified by finding a maximum matching in the bipartite graph. SCC consists of a spanning subgraph with the same node set as G, a maximum matching of G and an additional link set. The matching links divide the network G into several elementary paths and elementary cycles. The additional links transmit signals from the elementary paths to the elementary cycles. The unmatched nodes (in red) comprise the minimal set of driver nodes that can control the dynamics of the entire network from any initial state to any desired final state. b, c Overview of OptiCon for identifying optimal control nodes (OCNs) and their synergistic combinations. b Using gene expression data under two conditions (e.g., diseased vs. healthy) and a directed gene regulatory network as inputs to OptiCon, a deregulation score (DScore)-weighted network can be obtained. The control region of a gene in the DScore-weighted network consists of a direct control region and an indirect control region. Direct control region (highlighted in cyan) of a gene is identified by finding the structural control configuration of the network. Indirect control region (highlighted in yellow) is identified by using the indirect control value (ICV) and a shortest path (SP) search procedure. The candidate OCNs for combination therapy can be identified using a combinatorial optimization procedure. For clarity, only the control regions of g1, g5 and g7 are shown instead of all genes. o, d, and u, the optimal influence, desired influence and undesired influence by the candidate OCNs, respectively. c Identification of synergistic OCN pairs using synergy score. The synergy score consists of two parts. The mutation score measures the enrichment of recurrently mutated cancer genes in the optimal control region (OCR) of each OCN. The crosstalk score measures the interaction density between genes in the OCRs of the two OCNs under consideration. Norm, min-max normalization

Fig. 2

Performance assessment of predicted regulators. a Number of key regulators predicted by OptiCon, TargetControl, VIPER, RACS and a degree-based method. HCC, hepatocellular carcinoma; LUAD, lung adenocarcinoma; BRCA, breast invasive carcinoma. b Enrichment p-values of known cancer drug targets documented in the Therapeutic Target Database. RACS is not included in the comparison because RACS predictions are based on known cancer drug targets. c Enrichment p-values of side effect-causing proteins. No overlap with OptiCon predictions in HCC and LUAD. d Distribution of CERES scores of identified key regulators. P-values in d were computed using one-sided Kolmogorov–Smirnov test. The rest of the p-values were computed using hypergeometric distribution. Yellow dashed line, enrichment p-value of 0.05

Fig. 3

Performance assessment of predicted synergistic interactions. Top ranked synergistic gene pairs, predicted by OptiCon, VIPER, and RACS are evaluated using synergy score and enrichment of synthetic lethal interactions. a–c Synergy scores of predicted gene pairs. d–f Enrichment p-values of experimentally derived synthetic lethal interactions and g–i clinically relevant synthetic lethal interactions between the subnetworks targeted by predicted gene pairs. p-values for comparing methods in all panels were computed using one-sided Wilcoxon test. Enrichment p-values in panels (d–i) were computed using hypergeometric distribution. Lower, middle and upper lines of boxplots represent first quartile, median, and third quartile respectively. Lower and upper whiskers represent smallest values within 1.5 times interquartile range below Q1 and above Q3, respectively. Yellow dashed line, enrichment p-value of 0.05. HCC hepatocellular carcinoma, LUAD lung adenocarcinoma, BRCA breast invasive carcinoma. No synergistic pairs were predicted by RACS for HCC

Fig. 4

Synergistic optimal control nodes discovered in three cancer types. a, d, g 77, 192, and 63 significantly synergistic optimal control node (OCN) pairs (Benjamini–Hochberg adjusted empirical p-values < 0.05) identified in hepatocellular carcinoma (HCC), lung adenocarcinoma (LUAD) and breast invasive carcinoma (BRCA), respectively. Shade of red in the heat map is proportional to the synergy score. Numbers represent the ranks of the identified synergistic pairs. b, e, h Experimentally derived cancer-type-specific synthetic lethal interactions are significantly enriched between optimal control regions (OCRs) of 53 (69%), 138 (72%), and 61 (97%) synergistic OCN pairs identified in HCC, LUAD and BRCA, respectively (Benjamini–Hochberg adjusted hypergeometric p-values < 0.05). Shade of blue in the heat map is inversely proportional to the enrichment p-values. c, f, i Contingency tables and corresponding Fisher’s exact test p-values are shown for HCC, LUAD and BRCA, respectively, indicating that known synthetic lethal interactions (SL) are more enriched between OCRs of synergistic OCN pairs than non-synergistic OCN pairs

Fig. 5

Representative synergistic OCN pairs predicted for three cancer types. a Representative synergistic optimal control node (OCN) pair predicted for liver cancer, NCSTN and OGT, and crosstalk links between their specific optimal control regions (OCRs). For clarity, only genes that are known to be recurrently mutated in liver cancer and/or involved in known synthetic lethal (SL) interactions are labelled. Shade of a node represents the deregulation score (DScore) of the corresponding gene. Red, up-regulation in cancer samples; green, down-regulation in cancer samples. Yellow links, crosstalk links between the OCRs of the two OCNs. Blue dashed lines, known cancer-type-specific synthetic lethal interactions. b Representative synergistic OCN pair predicted for lung cancer, PARP1 and HIF1A, and crosstalk links between their specific OCRs. c Representative synergistic OCN pair predicted for breast cancer, PLK1 and PTP4A1, and crosstalk links between their specific OCRs

Fig. 6

Experimental validation of predicted synergistic pairs. a Representative FACS plot of cells infected with lentiviruses expressing OGT-targeting (GFP) and NCSTN-targeting (mCherry) sgRNAs. Day 0 (4 days post-transduction) and Day 5 data are shown. Value in each quadrant indicates the percentage of cells expressing a given reporter in the culture. The growth phenotype is calculated by measuring the relative depletion of the single-infected and double-infected cells between the start and the end of the growth assay. KO, single knockout. DKO, double knockout. RFU, relative fluorescence unit. Synergistic optimal control nodes validated in liver cancer (b), lung cancer (c), and breast cancer (d). Safe indicates non-targeting control sgRNA. GI score, genetic interaction score. Data represent mean ± s.d. from three replicate cultures. P-values were computed using one-sided t-test. Source data are provided as a Source Data file

Fig. 7

Crosstalk genes play an important role in therapy resistance. Effect of crosstalk genes on the interaction density between optimal control regions (OCRs) of an optimal control node (OCN) pair is quantified as the decrease in interaction density of two OCRs after a crosstalk gene is removed from the OCR, herein termed ΔD. For each cancer type, crosstalk genes with significant ΔD (empirical p-value < 0.1) are shown. Empirical p-value was calculated using a null distribution of crosstalk genes controlled by one million randomly selected gene pairs from the input network. Gene symbols are ordered from top to bottom in ascending statistical significance. 21 (50%), 26 (41%), and 10 (91%) of the crosstalk genes controlled by synergistic OCN pairs (magenta dots) have a known role in drug resistance in liver cancer (a), lung cancer (b) and breast cancer (c), respectively (Supplementary Data File 6). In contrast, 0, 7 (28%) and 2 (17%) of the crosstalk genes controlled by non-synergistic OCN pairs (yellow triangles) have a known role in drug resistance in the respective cancer types. HCC hepatocellular carcinoma, LUAD lung adenocarcinoma, BRCA breast invasive carcinoma