Control of Cellular Differentiation Trajectories for Cancer Reversion

Abstract

Cellular differentiation is controlled by intricate layers of gene regulation, involving the modulation of gene expression by various transcriptional regulators. Due to the complexity of gene regulation, identifying master regulators across the differentiation trajectory has been a longstanding challenge. To tackle this problem, a computational framework, single-cell Boolean network inference and control (BENEIN), is presented. Applying BENEIN to human large intestinal single-cell transcriptome data, MYB, HDAC2, and FOXA2 are identified as the master regulators whose inhibition induces enterocyte differentiation. It is found that simultaneous knockdown of these master regulators can revert colorectal cancer cells into normal-like enterocytes by synergistically inducing differentiation and suppressing malignancy, which is validated by in vitro and in vivo experiments.

Keywords: Boolean gene regulatory network model; cancer reversion; cell fate control; network reconstruction; single‐cell transcriptome; systems biology.

Figure 1

Schematic diagram of BENEIN for reconstruction of Boolean GRN from single‐cell transcriptome data and identification of master regulators for cancer reversion. A) BENEIN quantifies the abundance of pre‐mature and mature mRNA reads to separate the transcriptional status of each single‐cell into pre‐ and post‐transition states. B) BENEIN infers a potential regulatory structure between TFs and their TGs with a moving window strategy, by computing CMI and then using cisTarget database for eliminating indirect interactions for each window. C) BENEIN converts gene expressions of the pre‐ and post‐transition states into the binarized forms, where a value of 1 indicates that the gene is switched on and a value of 0 indicates that the gene is switched off. D–F) BENEIN infers Boolean functions for each gene by integrating the binarized matrices with the regulatory network structure to construct truth tables, then employing the QM algorithm. By iterating the process, BENEIN reconstructs the Boolean GRN model. G) BENEIN reduces the Boolean GRN model to its essential structure and finds minimal feedback vertex sets of the reduced network model to identify control target candidates. H) The identified control targets are optimized and validated by attractor simulation.

Figure 2

Inference and analysis of the Boolean GRN model for enterocyte differentiation. A) Single‐cell transcriptome data of enterocyte differentiation upon the UMAP embedding space with velocity stream plot (top) and pseudotime (bottom). B) The Boolean GRN model for enterocyte differentiation constructed by BENEIN. Red and blue edges represent activation and repression, respectively. C) Heatmap illustrating the gene expression dynamics along with the pseudotime. Gene expressions were binarized based on RNA splicing dynamics calculated by scVelo, then smoothed and converted to z‐score for improved visualization. Based on the gene expression dynamics, the binarized states for stem cell and enterocyte were determined by comparing the average gene activation in the cluster of each cell type and the average gene activation of all the cells in the trajectory. D) State transitions of the Boolean GRN model from a stem cell initial state (red) to an enterocyte attractor (blue) mapped onto the UMAP embedding space (see the Experimental Section for details on determining the position of each network state upon the UMAP embedding space). E) 32 point attractors of the Boolean GRN model mapped onto the UMAP embedding space, with the relative basin of attraction and Hamming distance between each attractor and binarized expression of the nearest cell. Similarity scores of the attractors with respect to F) the stem cell state and G) the enterocyte state upon the UMAP embedding space. The similarity score is defined by one minus the normalized Hamming distance between an attractor and the binarized state of the cell type.

Figure 3

Identification of optimal control targets for enterocyte differentiation and numerical simulations for the effects of control targets on the Boolean GRN model. A) Identification of control target candidates. The original network is shown with red edges representing activation and blue edges representing inhibition (left). The reduced network was obtained by applying the BNSimpleReduction algorithm to the original network (middle). The minimal FVS consisting of 5 nodes (HDAC2, MYB, SPDEF, PRDM1, and FOXA2) is marked in yellow in the reduced network (right). B) Optimization of the control targets. The bar plot shows cosine similarity between the desired attractor and the average activity of each perturbed network. The simultaneous inhibition of MYB, HDAC2, and FOXA2 yielded the highest similarity (≈0.97) and were chosen as the optimal control targets. C) Attractor and phenotype landscape of the unperturbed and each perturbed Boolean GRN model. All possible combinations within the optimal targets are considered. Numbers on the attractor landscape indicate indices of the attractors, and the area indicates a basin of attraction for each attractor. Red/blue in the phenotype landscape indicates a basin of attraction of the undifferentiated/differentiated enterocyte state, respectively. The perturbation corresponding to simultaneous inhibition of optimal control targets resulted in a 100% basin of attraction for enterocyte phenotype states. D,E) State transition trajectories between unperturbed and controlled networks. The graph on the top illustrates all possible state transitions, with the most probable trajectory highlighted in green, along with the transition probability. The graph below illustrates the most probable trajectory with updated node states according to the state transition. Color coding of the state ranges from red (stem cell state) to blue (enterocyte state). Without perturbation, the colon stem cell state can transition toward an undifferentiated state through 13 possible trajectories (D). When the control targets were regulated, the colon stem cell state would transition toward the enterocyte state through 2412 possible trajectories (E).

Figure 4

The underlying dynamical regulation mechanism of enterocyte differentiation induced by regulating the control targets of the Boolean GRN model. A) The transformation of the original network into a reduced network and its further rearrangement into a hierarchical network. The original network is shown with red edges representing activation and blue edges representing repression (top left). The reduced network was obtained by applying the BNSimpleReduction algorithm (top right). The reduced network is rearranged to visualize its hierarchical structure. Control targets are colored green in each of the networks, and their canalizing effects on the downstream nodes are highlighted in yellow on the edges (bottom). B) The Boolean GRN model visualized with the effectiveness (edges) and the weighted outdegree centrality (nodes). Color coding of each edge represents the effectiveness. Size of the node represents the weighted outdegree centrality, computed using effectiveness as weights (left). Bar plot showing the weighted outdegree of the effectiveness of each node, in descending order (right). The control targets show the highest weighted outdegree. C) Bar plot showing the synergy score for every possible three‐node perturbation. The nodes are each symbolized with a unique number and color. The combination of optimal control targets shows the highest synergy score. D) Comparison of the canalizing effect from the simultaneous perturbation of MYB, HDAC2, and FOXA2 in the Boolean GRN model with the in vitro transcript quantification on NCM‐460. The network is a subnetwork of the reduced network from Figure 4A, consisting of the canalized nodes. Nodes are colored red (activation) or blue (inhibition) according to the canalization effect. The in vitro transcript quantification results are shown in bar charts, with a scramble knockdown sample colored in gray and triple knockdown samples colored in red (up‐regulated) or blue (down‐regulated). Data are presented as the mean  ± SEM; n = 3 measurements (two‐tailed t‐test: *p < 0.05, **p < 0.01, ***p < 0.001) (NCM‐460 with scramble knockdown, C; NCM‐460 with simultaneous knockdown of MYB, HDAC2, and FOXA2, MHF).

Figure 5

The simultaneous perturbation of control targets inhibited proliferation of three colorectal cancer cells in vitro and in vivo. A) The growth curves of colorectal cancer cells (HCT‐116, CACO‐2, and HT‐29) after knockdown of control targets (top). Cell growth rate was analyzed by IncuCyte. Representative images of crystal violet staining of cells (bottom). Data are presented as the mean ± SEM; n = 3 replicates (two tailed t‐test: **p < 0.01; ***p < 0.001). B–D) HCT‐116, CACO‐2, and HT‐29 were injected to female athymic nude mice (Foxn1^nu/nu) and proliferation of cancer cells was observed in tumor‐bearing mice. Changes in the volume of three colorectal tumors 23 days after tumor injection (B). Photographs of tumors resected after sacrifice on day 23 (C). Tumor weight of a resected tumor (D). Data are presented as the mean ± SEM; n = 5 measurements (two tailed t‐test: ***p < 0.001).

Figure 6

Suppression of control targets reverts three colorectal cancer cells into normal‐like enterocytes by inactivating MYC and WNT pathways. A) Protein and mRNA expression levels of MYB, HDAC2, and FOXA2 in HCT‐116, CACO‐2, and HT‐29 cells (Three colorectal cancer cells with scramble knockdown, C; three colorectal cancer cells with simultaneous MYB, HDAC2, and FOXA2 knockdown, KD). B) Scatter plot on the upper layer shows the transcriptome data from three colorectal cancer cells and their respective reverted cells in three colors (Green, HCT‐116; Orange, CACO‐2; Purple, HT‐29). Scatter plot on the lower layer shows the transcriptome data from TCGA colorectal cancer (red) and their adjacent normal (blue), along with the gradient color of the background representing a support vector machine score for TCGA colon cancer expression. The x‐ and y‐axes represent normal colon cell and embryonic stem cell signature scores, respectively. C) Enrichment plot shows that the normalized enrichment scores (NESs) of normal enterocyte signature were increased by simultaneous knockdown of the control targets for three colon cancer cells (top). Box plot displays results of gene ontology analysis using genes up‐regulated by the simultaneous knockdown of the control targets. The up‐regulated genes are associated with the normal colon signatures (bottom). D) Protein abundances were monitored by western blotting analysis of the representative genes of colonic enterocytes (KRT19, KRT20, and VDR). E) Enrichment plots illustrate that the NESs of MYC (top) and WNT (bottom) signatures were decreased by the simultaneous knockdown of the control targets for three colon cancer cells. F) Protein abundances were monitored by western blotting analysis of the representative genes of MYC and WNT pathways (TCF1, MYC, and β‐catenin). Transcript quantification by qRT‐PCR is presented relative to that in the scramble shRNA condition. GAPDH was used as a loading control. Data are presented as the mean ± SEM; n ≥ 3 measurements (two‐tailed t‐test: *p < 0.05, **p < 0.01, ***p < 0.001).