Modeling differentiation-state transitions linked to therapeutic escape in triple-negative breast cancer

Abstract

Drug resistance in breast cancer cell populations has been shown to arise through phenotypic transition of cancer cells to a drug-tolerant state, for example through epithelial-to-mesenchymal transition or transition to a cancer stem cell state. However, many breast tumors are a heterogeneous mixture of cell types with numerous epigenetic states in addition to stem-like and mesenchymal phenotypes, and the dynamic behavior of this heterogeneous mixture in response to drug treatment is not well-understood. Recently, we showed that plasticity between differentiation states, as identified with intracellular markers such as cytokeratins, is linked to resistance to specific targeted therapeutics. Understanding the dynamics of differentiation-state transitions in this context could facilitate the development of more effective treatments for cancers that exhibit phenotypic heterogeneity and plasticity. In this work, we develop computational models of a drug-treated, phenotypically heterogeneous triple-negative breast cancer (TNBC) cell line to elucidate the feasibility of differentiation-state transition as a mechanism for therapeutic escape in this tumor subtype. Specifically, we use modeling to predict the changes in differentiation-state transitions that underlie specific therapy-induced changes in differentiation-state marker expression that we recently observed in the HCC1143 cell line. We report several statistically significant therapy-induced changes in transition rates between basal, luminal, mesenchymal, and non-basal/non-luminal/non-mesenchymal differentiation states in HCC1143 cell populations. Moreover, we validate model predictions on cell division and cell death empirically, and we test our models on an independent data set. Overall, we demonstrate that changes in differentiation-state transition rates induced by targeted therapy can provoke distinct differentiation-state aggregations of drug-resistant cells, which may be fundamental to the design of improved therapeutic regimens for cancers with phenotypic heterogeneity.

Fig 1. Drug-specific model.

Live cells occupy four differentiation states and can transition, divide, or die. The dynamics parameters {ρ_ij, ρ_i, ρ_iD}_i≠j are the average rates of these actions taken by live cells in each differentiation state following treatment.

Fig 2. Ensemble model predictions in comparison to training data.

The training samples (black stars) and predictions by the model ensemble (gray bands) are shown for each treatment condition: DMSO (row 1), Trametinib (row 2), BEZ235 (row 3), and Trametinib+BEZ235 (row 4). The model ensemble is a collection of models that were identified from the training data via resampling residuals bootstrap [30] for each treatment condition. In each plot, we show a 95% confidence interval (gray band) around the median (black dotted line) of the ensemble model predictions. Higher p-values indicate better consistency between predictions and training data over the time horizon (12h, 24h, …, 72h).

Fig 3. Drug-specific transition gains.

For each treatment condition, values of the transition gains from the AM-optimized dynamics matrix are shown (units: $\frac{\#\text{cells}\text{at}(k+1)\text{multiples}\text{of}12\text{hours}}{\#\text{cells}\text{at}k\text{multiples}\text{of}12\text{hours}}$). Each transition gain from differentiation state i to differentiation state j of sufficient magnitude (ρ_ij ≥ 0.10) is depicted as an arrow directed from i to j. Arrow style specifies gain magnitude. A dotted arrow means ρ_ij ∈ [0.10, 0.30), a dashed arrow means ρ_ij ∈ [0.30, 0.70), and a solid arrow means ρ_ij ∈ [0.70, 1.00].

Fig 4. Uncertainty analysis of the dynamics parameters.

For each treatment condition, 95% confidence intervals computed from the model ensemble are shown. These intervals indicate variations of the dynamics parameters due to measurement noise. Non-overlapping intervals of a given parameter specify a statistically significant difference. For example, a statistically significant reduction in K14^hi-to-VIM^hiK14^low transition was detected under Trametinib versus DMSO because the ρ₁₂-interval for Trametinib is strictly below the ρ₁₂-interval for DMSO. A p-value for each dynamics parameter is also provided in S2 Appendix.

Fig 5. Further investigations of Trametinib-induced K14^hi enrichment hypothesis and Trametinib+BEZ235-induced VIM^hiK14^low de-enrichment hypothesis.

Top left: Trametinib K14^hi live cell predictions by the AM-optimized dynamics matrix (pink band) are shown in comparison to test data (black stars). Top right: Trametinib K14^hi live cell predictions by a dynamics matrix identified with additional constraints (pink band) in comparison to test data (black stars). The additional constraints are ρ₁₂ ≥ 0.59, the value of ρ₁₂ for DMSO, and ρ₃₁ ≤ 0.19, the value of ρ₃₁ for DMSO. (ρ₁₂ is the K14^hi-to-VIM^hiK14^low transition gain, and ρ₃₁ is the K19^hiVIM^lowK14^low-to-K14^hi transition gain). Bottom left: Trametinib+BEZ235 VIM^hiK14^low live cell predictions by the AM-optimized dynamics matrix (pink band) are shown in comparison to test data (black stars). Bottom right: Trametinib+BEZ235 VIM^hiK14^low live cell predictions by a dynamics matrix identified with additional constraints (pink band) are shown in comparison to test data (black stars). In this setting, the additional constraints are ρ₁₂ ≥ 0.59, the value of ρ₁₂ for DMSO, and ρ₃₂ ≥ 0.63, the value of ρ₃₂ for DMSO. (ρ₁₂ is the K14^hi-to-VIM^hiK14^low transition gain, and ρ₃₂ is the K19^hiVIM^lowK14^low-to-VIM^hiK14^low transition gain). In each plot, the pink band extends between the maximum prediction and the minimum prediction out of four predictions in total at each time point (0h, 12h, …, 72h). The dotted line indicates the median of the predictions. Higher p-values indicate better consistency between predictions and test data over the time horizon (12h, 24h, …, 60h).

Fig 6. Death time series data for model training.

The YO-PRO-1 dye was used to quantify the proportion of dying cells every 12h in response to drug treatment [18]. The cells were treated with DMSO (baseline), 1μM Trametinib, 1μM BEZ235, or the combination of 1μM Trametinib + 1μM BEZ235.

Fig 7. Differentiation-state time series portion of the test data.

The sample mean and the sample standard deviation of fold change for each differentiation state are shown at each time point, computed from 4 samples. Fold change is $\frac{\text{fraction}\text{differentiation}\text{state}i,\text{well}w,\text{time}k,\text{therapy}}{\text{average}\text{fraction}\text{differentiation}\text{state}i,\text{time}k,\text{DMSO}}$, where fraction differentiation state i is the number of cells counted in that state divided by the population total.

Fig 8. Ensemble model predictions in comparison to test data.

The test samples (black stars) and ensemble model predictions (gray bands) are shown for each treatment condition: DMSO (row 1), Trametinib (row 2), BEZ235 (row 3), and Trametinib+BEZ235 (row 4). The model ensemble is a collection of models that were identified from the training data via resampling residuals bootstrap [30] for each treatment condition. In each plot, we show a 95% confidence interval (gray band) around the median (black dotted line) of the ensemble model predictions. Higher p-values indicate better consistency between predictions and test data over the time horizon (12h, 24h, …, 60h).

Fig 9. Single model predictions in comparison to test data, where the differentiation states are defined by K14 only.

The test samples (black stars) and single model predictions (pink bands) are shown for each treatment condition: DMSO (row 1), Trametinib (row 2), BEZ235 (row 3), and Trametinib+BEZ235 (row 4). The single model was identified on the training data using K14^hi and K14^low as the differentiation states for each treatment condition. The pink band extends between the maximum prediction and the minimum prediction out of four predictions in total at each time point (0h, 12h, …, 72h). The dotted line indicates the median of the predictions. Higher p-values indicate better consistency between predictions and test data over the time horizon (12h, 24h, …, 60h).