Skip to main page content
U.S. flag

An official website of the United States government

Dot gov

The .gov means it’s official.
Federal government websites often end in .gov or .mil. Before sharing sensitive information, make sure you’re on a federal government site.

Https

The site is secure.
The https:// ensures that you are connecting to the official website and that any information you provide is encrypted and transmitted securely.

Access keys NCBI Homepage MyNCBI Homepage Main Content Main Navigation
. 2021 Sep 16;11(1):18511.
doi: 10.1038/s41598-021-97887-z.

A combination approach of pseudotime analysis and mathematical modeling for understanding drug-resistant mechanisms

Affiliations

A combination approach of pseudotime analysis and mathematical modeling for understanding drug-resistant mechanisms

Shigeyuki Magi et al. Sci Rep. .

Abstract

Cancer cells acquire drug resistance through the following stages: nonresistant, pre-resistant, and resistant. Although the molecular mechanism of drug resistance is well investigated, the process of drug resistance acquisition remains largely unknown. Here we elucidate the molecular mechanisms underlying the process of drug resistance acquisition by sequential analysis of gene expression patterns in tamoxifen-treated breast cancer cells. Single-cell RNA-sequencing indicates that tamoxifen-resistant cells can be subgrouped into two, one showing altered gene expression related to metabolic regulation and another showing high expression levels of adhesion-related molecules and histone-modifying enzymes. Pseudotime analysis showed a cell transition trajectory to the two resistant subgroups that stem from a shared pre-resistant state. An ordinary differential equation model based on the trajectory fitted well with the experimental results of cell growth. Based on the established model, it was predicted and experimentally validated that inhibition of transition to both resistant subtypes would prevent the appearance of tamoxifen resistance.

PubMed Disclaimer

Conflict of interest statement

The authors declare no competing interests.

Figures

Figure 1
Figure 1
Functional analysis of gene expression patterns in human breast adenocarcinoma MCF-7 cells during the TAM resistance acquisition process. (a) Schematic overview of the experimental procedure. Each treatment was replicated five times. (b) Growth rate of TAM-treated and control (Ctrl) cells. Data represent mean ± standard error (SE, n = 3). (c,d) Ratio of S (c) and G1 (d) phase in TAM-treated and Ctrl cells. Data represent mean ± SE (n = 3 in week 1 to week 6, n = 2 in week 7 to week10). From (b) to (d), the data points are slightly shifted in the x-axis direction to improve the visibility and prevent overlapping of the Ctrl and TAM graphs, and p-values were calculated using two-tailed Welch's test. (e,f) Cluster analysis of the z-scores of log2 fold change (log2FC) values by time points (e) and by genes (f). The bottom line graphs in (f) showed individual (orange) or median (red) of gene expression patterns. (g,h) Heatmaps of enrichment analysis data. The top five significant terms in the Reactome pathway database (g) and KEGG pathway database (h) in each cluster are presented.
Figure 2
Figure 2
Single-cell RNA-seq analysis of TAM-resistant MCF-7 cells. (a) Schematic overview of the experimental procedure of single-cell RNA-seq. (b) Boxplot showing the distribution of the correlation coefficients of single-cell gene counting among cells at each time point. The median values are presented in red. The q-values were calculated using Wilcoxon rank-sum test. (c) Percentage of cells at different cell cycle stages at each time point. (df) Visualization of single-cell transcriptome data by UMAP. Single-cell data space was reduced to three dimensions, and the distribution of data was visualized using the first two dimensions. Cells were colored by estimated cell cycle stage (d), week (e), and subgroups (f). (g) Complex heatmap of cell subpopulation, co-regulated gene modules, and enriched functions. Top left heatmap presents the frequency of subgroups at each time point shown in (f). Top right heatmap presents the relative expression level of each gene module in each cluster. Bottom heatmap presents the enrichment terms in each gene module.
Figure 3
Figure 3
Trajectory analysis of the TAM resistance acquisition process. (ac) Single-cell trajectory during the continuous TAM treatment. Graphs were colored by clusters (a), weeks (b), and pseudotime (c). Numbered circles in white and black indicate a root node and estimated branch nodes in the trajectory, respectively. (d,e) Violin plots showing the distribution of pseudotime in each week (d) and subgroup (e). The q-values in (d) are calculated using Wilcoxon rank-sum test. (f) Expression score of up-regulated genes in subgroup 4 (left) and subgroup 5 (right), considering the cell trajectory. (g) Upstream factor analysis of up-regulated genes displayed in (f). Top 10 q-value data are presented. Color represents subtypes of breast cancers: yellow, MCF-7 cells; green, triple-negative breast cancer (TNBC). (h) Time-series bulk gene expression patterns listed in (g). NELFE was not detected in bulk RNA-seq.
Figure 4
Figure 4
Ordinary differential equation-based model of the TAM resistance acquisition process. (a) Illustration of the model scheme. (b,c) Time-series analysis of changes in cell growth rate (b) and ratio of each subpopulation (c). Points: experimental data; error bars: standard deviation (SD) of experimental data; solid lines: averaged in silico simulation of 20 sets of parameters; shaded areas: SD of simulations. (d) Results of sensitivity analysis of the mean growth rate from week 3 (W3) to week 10 (W10) at each reaction; v1v12. Error bars represent the SD of simulations with 20 set parameters. The p-value was calculated using Wilcoxon signed-rank test. (e) Heatmap of the mean growth rate of TAM-treated cells from W3 to W10, with different inhibitory intensity of growth rate of XR1 or XR2 and cell transition to XR1 or XR2. The X- and Y-axes indicate the remaining reaction rate (i.e., 1.0 means 0% inhibition, and 0.0 means 100% inhibition). Data represent the mean of simulations with 20 set parameters. Blue lines represent a border dividing mean growth rate of less than or equal to 1. Gray boxes indicate where conditions were compared in the main text and show q-values < 0.05.
Figure 5
Figure 5
Experimental validation. (a) Effect of PML knockdown and KDM5 inhibitor GSK467 on the relative cell number of TAM-treated cells measured using MTT assay at each time point. siCtrl and siRNA indicate conditions in non-targeting siRNA-treated and those in siRNA targeting PML-treated, respectively. The data represent mean ± SE (n = 4). q-values were calculated using two-tailed Welch's test. (b) Protein expression levels of PML in each condition. β-actin presented as a loading control. Whole membrane images are presented in Supplementary Fig. S9.

Similar articles

Cited by

References

    1. Aranda A, Pascual A. Nuclear hormone receptors and gene expression. Physiol. Rev. 2001;81:1269–1304. doi: 10.1152/physrev.2001.81.3.1269. - DOI - PubMed
    1. Burns KA, Korach KS. Estrogen receptors and human disease: An update. Arch. Toxicol. 2012;86:1491–1504. doi: 10.1007/s00204-012-0868-5. - DOI - PMC - PubMed
    1. Allred DC, Brown P, Medina D. The origins of estrogen receptor alpha-positive and estrogen receptor alpha-negative human breast cancer. Breast Cancer Res. 2004;6:240–245. doi: 10.1186/bcr938. - DOI - PMC - PubMed
    1. Ring A, Dowsett M. Mechanisms of tamoxifen resistance. Endocr. Relat. Cancer. 2004;11:643–658. doi: 10.1677/erc.1.00776. - DOI - PubMed
    1. Jeselsohn R, Buchwalter G, De Angelis C, Brown M, Schiff R. ESR1 mutations—a mechanism for acquired endocrine resistance in breast cancer. Nat. Rev. Clin. Oncol. 2015;12:573–583. doi: 10.1038/nrclinonc.2015.117. - DOI - PMC - PubMed

Publication types

MeSH terms