Modeling glioblastoma heterogeneity as a dynamic network of cell states

Abstract

Tumor cell heterogeneity is a crucial characteristic of malignant brain tumors and underpins phenomena such as therapy resistance and tumor recurrence. Advances in single-cell analysis have enabled the delineation of distinct cellular states of brain tumor cells, but the time-dependent changes in such states remain poorly understood. Here, we construct quantitative models of the time-dependent transcriptional variation of patient-derived glioblastoma (GBM) cells. We build the models by sampling and profiling barcoded GBM cells and their progeny over the course of 3 weeks and by fitting a mathematical model to estimate changes in GBM cell states and their growth rates. Our model suggests a hierarchical yet plastic organization of GBM, where the rates and patterns of cell state switching are partly patient-specific. Therapeutic interventions produce complex dynamic effects, including inhibition of specific states and altered differentiation. Our method provides a general strategy to uncover time-dependent changes in cancer cells and offers a way to evaluate and predict how therapy affects cell state composition.

Keywords: cell state; cellular barcoding; patient-derived brain tumor cells; single-cell lineage tracing; time-dependent computational models.

Figure 1. Resolving the plasticity of GBM cells

1. A

   Overview. In the STAG procedure, barcoded tumor cells are sampled at multiple time points and profiled by single‐cell RNA sequencing, followed by mathematical modeling to identify (i) the states and their growth rates, (ii) the patterns and rates of the state transitions, (iii) how drugs affect cell states, and (iv) analysis of cell population stability and long‐term projections of cell state compositions.

2. B

   Clone sizes of barcoded U3065MG cells at 0–21 days. X‐axis, barcodes; Y‐axis, number of cells containing each barcode.

3. C

   Venn diagram of number of barcodes detected at 7, 14, and 21 days.

4. D, E

   Simulation of the experiment in (B,C) using a Galton–Watson process with a fixed growth and death rate.

5. F

   Enriched gene sets in small and large clones, respectively, using Gene Set Enrichment Analysis (GSEA). Normalized enrichment scores (NES) and q‐values are indicated in the figure.

Figure 2. State dynamics of cells from the mesenchymal U3065MG GBM cell line

1. UMAP embedding of single U3065MG cells colored according to transcriptional cell state assignment.

2. Average number of states represented in each clone with increasing clone size. X‐axis, clone size; Y‐axis, average number of states.

3. UMAP embedding of cell state distribution at each experimental time point. State colors as in (A).

4. Heatmap of the U3065MG state transition and growth matrix A derived from the STAG model.

5. Estimated network of state transitions in the U3065MG cell line. The thickness of the arrows correlates with the rate of transition (number of transitions per day).

6. STAG estimates of growth rates (cell divisions per day) for U3065MG states 1–6. Error bars indicate the 90% confidence interval, based on 1,000 bootstrap replicates.

7. Enrichment of the U3065MG state 1 cell population (CD24^high/CD44^high/HLA‐DR^low) using FACS cell sorting. Cells were monitored by flow cytometry over 3 weeks, and a gradual phenotypic drift, e.g., toward state 4 (CD24^high/CD44^high/HLA‐DR^high), could be observed. See also Fig EV1.

8. scRNA sequencing of a U3065MG‐derived clonal culture (clone 475) detected all six states of the U3065MG mother cell line. U3065MG‐C475 cells were scored against U3065MG state 1–6 signatures using ssGSEA.

Figure EV1. Cell surface marker selection and extended FACS analyses

1. Heatmap of selected cell surface marker mRNA expression (Z scores) in each of the U3065MG states.

2. Flow cytometry analysis of cell surface marker expression, 2 weeks after FACS‐enrichment of U3065MG state 1 cells (CD24^high/CD44^high/HLA‐DR^low), as in Fig 2G. Two biological replicate experiments are shown.

Figure 3. Characterization of states in the mesenchymal U3065MG GBM cell line

1. GSEA of the gene signatures for each of the six U3065MG states. Blue and red indicate negative and positive enrichments of the pathway. Size of the dot indicates significance (proportional to the GSEA false discovery rate [FDR] q‐value). Gene sets without PMID reference were obtained from the MSigDB Hallmarks database of gene sets. Abbreviations: oligodendrocyte progenitor cell (OPC), neural progenitor cell (NPC), mesenchymal (MES), astrocyte (AC), neuronal (NEU), mitochondrial (MTC), proliferative/progenitor (PPR).

2. Survival analysis of TCGA GBM patients estimated by Cox's proportional hazards model, with enrichment score of states 1–6 as independent covariate. Shaded areas indicate 95% confidence intervals, calculated as ± 1.96 * standard error (SE)). HR, hazard ratio.

Source data are available online for this figure.

Figure EV2. Cluster robustness analysis U3065MG

Sankey plot showing how cells change state assignment when k increases from 4 to 8. Inspecting the clustering solution for k between 4 and 8, the main trend is that increasing k by 1 tends to split a cluster in two, indicating that clusters are robust.

Figure 4. State dynamics of cells from the mesenchymal U3071MG and the classical U3017MG GBM cell lines

1. A, B

   UMAP embedding of single cells from the mesenchymal U3071MG (A) and classical U3017MG (B) GBM cell lines, colored according to transcriptional state assignments I–VI and A–F, respectively.

2. C, D

   Estimated network of state transitions in U3071MG (C) and U3017MG (D). The thickness of the arrows correlates with the rate of transition (number of transitions per day, as in Fig 2E).

3. E‐F

   State growth rates (number of cell divisions per day) in U3071MG (E) and U3017MG (F).

4. G, H

   Sankey plot of the relation between the states defined in U3065MG cells (1–6), and the states defined in U3071MG cells (I–VI) (G) and in U3017MG cells (A–F) (H). U3071MG and U3017MG cells were scored against U3065MG state 1–6 signatures using ssGSEA.

Figure 5. Characterization of the states in the mesenchymal U3071MG and the classical U3017MG GBM cell lines

1. GSEA of the gene signatures for each of the six states in U3017MG and U3071MG. Blue and red indicate negative and positive enrichments of the pathway. Size of the dot indicates significance (proportional to the GSEA FDR q‐value). Gene sets without PMID reference were obtained from the MSigDB Hallmarks database of gene sets. Abbreviations: oligodendrocyte progenitor cell (OPC), neural progenitor (NPC), mesenchymal (MES), astrocyte (AC), neuronal (NEU), mitochondrial (MTC), proliferative/progenitor (PPR).

2. State signatures for U3065MG (1–6), U3017MG (A–F), and U3071MG (I–VI) related to embryonic cluster definitions from Eze et al. The heatmap shows enrichment results; red and blue indicate positive and negative enrichments, respectively. The dot plot in the right panel shows mean Carnegie stage for cells in each embryonic cluster. Barplots below the heatmap show Spearman correlation (n = 61) between the enrichment profile for each state and the Carnegie time for each embryonic cluster. Spearman P‐values are shown below barplots.

Source data are available online for this figure.

Figure EV3. Survival analysis U3017MG and U3071MG

State signatures for states in U3017MG and U3071MG were related to the TCGA GBM cohort through ssGSEA and enrichment scores were used as independent covariate to build Cox's proportional hazards models. Shaded areas indicate 95% confidence intervals, calculated as ± 1.96 * standard error (SE).

Figure 6. Employing the transition networks for assessing and predicting therapeutic interventions

1. Experimental design for U3065MG treatment experiments.

2. STAG estimates of the state‐specific growth rates for untreated U3065MG cells, and cells treated with BMP4 and TMZ. Error bars indicate the 90% confidence interval, based on 1,000 bootstrap replicates.

3. Estimated network of state transitions in the untreated U3065MG cells. The thickness of the arrows correlates with the rate of transition (number of transitions per day).

4. Estimated network of state transitions for BMP4‐ and TMZ‐treated U3065MG cells. Red indicates increased rate of transitions compared with untreated cells. Blue indicates decreased rate of transitions compared with untreated cells. Gray indicates transitions of the untreated network.

5. Conceptual image of how the STAG‐derived A‐matrix can be expressed as eigenvalues and eigenvectors, interpreted as measures of population stability and state equilibria, respectively.

6. Eigenvalues of the state transition network of TMZ‐treated U3065MG cells and the predicted minimal intervention to obtain bounded growth in each case.

7. Predicted state equilibria for U3065MG cells in untreated condition and treated with BMP4 or TMZ.

Figure EV4. State equilibria U3017MG and U3071MG

Stacked barplots showing the predicted state equilibria for U3017MG states A–E and U3071MG states I–VI at steady state. The state F‐fraction in U3017MG is < 0.001 and not shown in figure.

Figure EV5. Cluster robustness analysis U3017MG and U3071MG

Sankey plot showing how cells change state assignment when k increases from 4 to 8.