www.nature.

com/npjsba

ARTICLE OPEN

Understanding and leveraging phenotypic plasticity during

Cancer metastasis is the process of detrimental systemic spread and the primary cause of cancer-related fatalities. Successful
metastasis formation requires tumor cells to be proliferative and invasive; however, cells cannot be effective at both tasks
simultaneously. Tumor cells compensate for this trade-off by changing their phenotype during metastasis formation through
phenotypic plasticity. Given the changing selection pressures and competitive interactions that tumor cells face, it is poorly
understood how plasticity shapes the process of metastasis formation. Here, we develop an ecology-inspired mathematical model
with phenotypic plasticity and resource competition between phenotypes to address this knowledge gap. We find that
phenotypically plastic tumor cell populations attain a stable phenotype equilibrium that maintains tumor cell heterogeneity.
Considering treatment types inspired by chemo- and immunotherapy, we highlight that plasticity can protect tumors against
interventions. Turning this strength into a weakness, we corroborate current clinical practices to use plasticity as a target for
adjuvant therapy. We present a parsimonious view of tumor plasticity-driven metastasis that is quantitative and experimentally
testable, and thus potentially improving the mechanistic understanding of metastasis at the cell population level, and its treatment
1234567890():,;

consequences.
npj Systems Biology and Applications (2023)9:48 ; https://doi.org/10.1038/s41540-023-00309-1

INTRODUCTION cannot move, while mesenchymal cells are motile and invasive
Cancer metastasis formation, or the spread and growth of tumor but proliferate slowly12. Both proliferative and invasive pheno-
cells throughout the body, is the cause of more than 60% of types are crucial for tumor progression and metastasis. Thus, shifts
cancer-related deaths1. Despite many decades of drug develop- in the phenotype of tumor cells by epithelial-mesenchymal
ment for cancer, the survival rate for patients with metastatic plasticity are salient features of metastatic cancers13–16. Conse-
cancer remains low2,3. The process of metastasis is a multistage quently, during the epithelial-to-mesenchymal transition, partial
process that includes local invasion by the tumor cells, invading mesenchymal features are gained, and epithelial features are
into the circulatory system, survival in circulation, arrest at a progressively lost, leading to altered surface marker expression,
distant tissue, exiting the circulation, survival, adaptation and decreased cell-cell adhesion, and motility-facilitating cytoskeleton
outgrowth in a new environment4–6. Tumor cells face selection for reorganization17. The reverse change happens during the
proliferation locally at each site of cancer, and for invasiveness and mesenchymal-to-epithelial transition. These transitions result in
motility during the spread between two sites. These changing the presence of hybrid phenotypes along the continuous
selection pressures during metastasis are driven by environmental spectrum between epithelial and mesenchymal phenotypes. This
factors and interactions with different surroundings7,8. continuous spectrum is often discretized into a number of hybrid
How tumor cells respond to these frequently changing phenotypes, with varying estimates on the number of these
environmental conditions is crucial for their persistence. Tumor hybrid phenotypes18–22. The plasticity between epithelial and
cell populations possess several adaptive strategies to cope with mesenchymal phenotypes was found to be promoted by
such fluctuating environments. A subset of these strategies relies interactions with the tumor microenvironment. Among others,
on diversity-generating mechanisms. In addition to genetic high levels of Transforming Growth Factor β (TGF-β) is a critical
mechanisms, such as mutations, copy number alterations, and driving factor toward more mesenchymal phenotypes8,17,23.
translocations, there are non-genetic mechanisms to generate Different treatment types select differently on this plasticity24.
diversity, such as epigenetic regulation, stochastic gene expres- The current first-line treatment for many cancers is chemotherapy
sion, and cellular differentiation hierarchies7,9. All these factors are which targets cell proliferation. A more recent treatment option
internal to the cell. In addition, the effect of the local environment for some cancers is immunotherapy which comprises antibody-
can be substantial; local variations in the tumor microenvironment based or cellular therapies targeting molecules present on the cell
may generate functional heterogeneity among clonal tumor membrane of the cancer cells. These differential selection
cells7,8. An early review by West-Eberhard10 emphasizes plasticity pressures suggest classifying treatments into growth-dependent
as a diversity-generating mechanism and its contribution to or growth-independent treatments. As epithelial-like phenotypes
altering traits. Plasticity also affects population and eco- proliferate faster, they are more sensitive to chemotherapy than
evolutionary dynamics11. mesenchymal phenotypes. Immunotherapy should target all
One prominent outcome of epithelial tumor cell plasticity is the phenotypes similarly, assuming that it does not target surface
range of phenotypes and transitions on the spectrum from markers that change during epithelial-mesenchymal transition.
epithelium to mesenchyme. Epithelial cells are proliferative but Adjuvant therapies, such as TGF-β blockers, are promising

1
Department of Theoretical Biology, Max Planck Institute for Evolutionary Biology, August-Thienemann-Str. 2, 24306 Plön, Germany. 2Institute for Experimental Cancer Research,
Kiel University and University Hospital Schleswig-Holstein, Campus Kiel, Arnold-Heller-Str. 3, Building U30, Entrance 1, 24105 Kiel, Germany. ✉email: [email protected]

Published in partnership with the Systems Biology Institute

Hybrid

Mesenchymal
Epithelial

Spread
Primary site

Secondary site

Fig. 1 Model structure capturing phenotypic heterogeneity at the primary and secondary tumor site. We focus on the competitive growth
of heterogeneous tumor cell populations at each tumor site. The dashed arrow represents spread of cancer cells between the primary and
1234567890():,;

secondary sites. We do not model the spread between tumor sites explicitly, but consider the spread as translating into different initial
phenotype distributions at primary and secondary sites. The different compartments in the model for each site represent epithelial, hybrid,
and mesenchymal phenotypes. The solid arrows indicate competitive growth and phenotype transitions. Phenotype i grows at rate ri and
transitions to the adjacent more mesenchymal-like type at rate TEM and to the adjacent more epithelial-like type at rate TME. Resource
P
competition is modeled with the term r i XK where X ¼ Nj¼1 x j is the total population abundance, and K is the carrying capacity. The epithelial
phenotype can only transition to the more mesenchymal-like adjacent hybrid phenotype, and the mesenchymal phenotype can only
transition to the epithelial-like adjacent hybrid phenotype; thus, the epithelial and mesenchymal phenotypes are the terminal phenotypes.
Here, we show only one hybrid phenotype, but we also investigate the effect of a larger number of hybrid phenotypes.

candidates to modify the transition process and, thus, ultimately, heterogeneous response to drugs as an inherent property of the
the phenotypic plasticity of tumor cells25–27. As TGF-β promotes tumor. We also discuss feasible observables that could be
the mesenchymal phenotype26,27 and TGF-β blockers improve measured in patient-derived tumors to test our hypothesis. In
chemotherapy25, it is suspected that the gene regulatory network this study, we focus on (i) how different characteristics of plasticity
of the epithelial-mesenchymal transition could be a key to control shape the phenotype distribution, (ii) how this distribution
the balance of phenotypes and thus the treatment outcome28,29. determines the transition dynamics between epithelial and
This gene regulatory network consists of cell fate-controlling mesenchymal phenotypes at primary and secondary tumor sites,
signaling pathways like TGF-β, WNT, and NOTCH that mediate the and (iii) how these dynamics are affected by chemo- or
recruitment of epithelial-mesenchymal-transition-inducing tran- immunotherapy, and additionally by a transition-modulating
scription factors ZEB, SNAIL, SLUG, TWIST, or miR-20017,29–32.
therapy. With the proposal of a practical transition modulating
Mathematical modeling presents a valuable approach to
therapy we emphasize the phenotypic plasticity as a viable clinical
understand the mechanisms and consequences of genetic and
target.
non-genetic diversity-generating mechanisms. Along this line
Zhou et al.33, modeled the hierarchical cellular diversity created
by differentiation and de-differentiation, with stem cells at the RESULTS
root of the hierarchy and specialized cells at its leaves. Gupta Overview
et al.34 and Li and Thirumalai35 showed with experiments and
mathematical modeling that phenotypically plastic breast cancers Using a dynamical model of the abundances, xi, of N distinct
approach a phenotypic equilibrium and can maintain phenotypic phenotypes with linearly decreasing growth rates along the
heterogeneity. A mathematical model of phenotype switching epithelial-mesenchymal trait axis, i = 1, 2, …, N, we track the
between drug-sensitive and drug-resistant cells in non-small cell temporal phenotype dynamics before, during and after treatment
lung cancer36 also shows a similar behavior. Such plastic (see Fig. 1). We implement competition between phenotypes by
mechanisms in tumor cells cause reduced efficacy of chemother- assuming a shared carrying capacity K. We include transitions to
apy and acquired therapeutic resistance13,37,38. In a spatial, multi- an adjacent mesenchymal-like phenotype TEM and adjacent
organ model, Franssen and Chaplain39 focus specifically on epithelial-like phenotype TME. We propose that a factor can
epithelial-mesenchymal transitions during metastatic spread with- control the balance between these transitions TEM and TME. The
out treatment. In a recent paper, we modeled the phenotype concentration of this transition-modulating factor can be mapped
switching between fast and slowly proliferating B cells during the to the transition bias λ. The relation between growth and speed of
relapse dynamics in acute lymphoblastic leukemia and investi- phenotype transition is encoded by transition speed c. Further, we
gated the impacts of treatment on phenotypic heterogeneity24. investigate growth-dependent treatment mD, i.e., chemotherapy,
In the present study, we propose an unconform yet parsimo- and growth-independent treatment mI, i.e., immunotherapy, to
nious view of the phenotype transitions without explicitly explore treatments and administration sequences. For a detailed
considering the microenvironment. This view considers the trait description of the modeling assumptions please see the Section
variation present in a tumor population and consequently the Model.

npj Systems Biology and Applications (2023) 48 Published in partnership with the Systems Biology Institute
 S. Shah et al.
3

Fig. 2 Phenotypic plasticity creates a stable phenotype distribution. Each panel shows the distribution of phenotypes (Eq. (1)) relative to
the carrying capacity K for a fixed value of the transition bias λ and the number of phenotypes N. The stable phenotype distribution changes
with the transition bias λ but remains qualitatively unaffected by changing the number of phenotypes N. The stable distribution is uniform
when there is no transition bias to either epithelial or mesenchymal-like phenotypes, i.e., λ = 0. λ < 0 depicts a transition bias towards
epithelial-like phenotypes and leads to a relative increase in epithelial cells. Conversely, λ > 0 results in a transition bias towards mesenchymal-
like phenotypes and causes a relative increase in mesenchymal cells.

Phenotype transitions generate and maintain heterogeneity Transition bias controls the stable phenotype distribution
Assuming that cells can transition into adjacent epithelial or When cells change to a more epithelial or mesenchymal
mesenchymal phenotypes generates a phenotype distribution of phenotype with the same probability, i.e., when there is no
abundances along the epithelial-mesenchymal trait axis. Without transition bias (λ = 0), all phenotypes are at equal abundance and
treatment, mD = mI = 0, the model (see Section Model) has two the stable distribution is a uniform distribution (Fig. 2). The
equilibria. One equilibrium is the state where the tumor vanishes, mesenchymal phenotype (M) is the most abundant when there is
xi = 0 for i = 1, 2, …, N. This equilibrium is unstable, implying that a a transition bias towards the mesenchymal phenotypes (λ > 0). On
tumor can always progress and increase in abundance without the other hand, for a transition bias towards the epithelial
treatment. The second equilibrium is the state described by phenotypes (λ < 0), the epithelial phenotype (E) is the most
abundant. The phenotypic heterogeneity of the tumor population
ð1  λÞNi ð1 þ λÞi1 is highest when there is no transition bias and λ = 0 (Supplemen-
x i ¼ K PN Nj
for i ¼ 1; 2; ¼ ; N: (1) tary Fig. 1). A high variance corresponds to the presence of a lot of
j¼1 ð1  λÞ ð1 þ λÞj1 different phenotypes reducing the efficacy of standard cancer
treatments. The heterogeneity decreases and vanishes at the bias
This equilibrium is stable and describes the coexistence of all extremes λ → ± 1 where transitions become unidirectional.
phenotypes (see Supplementary Section Stability of the coex- The presence of phenotypes far off the average population
istence equilibrium for a proof of global stability). This stability phenotype with very low abundances is captured by the third
implies that phenotypic heterogeneity is generated and main- central moment of the phenotype distribution. These phenotypes
tained in tumor populations featuring phenotype transitions. are essential for residual disease in some cancers, and thus an
Patient-derived breast cancer cell populations exhibit similar important potential treatment target. We find that the third
behavior in-vitro34,35. central moment is an odd function of the transition bias, vanishing
The sum of all phenotype abundances at the coexistence at λ = 0 and at high transition bias (Supplementary Fig. 1). The
equilibrium is K. The stable distribution of phenotypes at the number of phenotypes N quantitatively affects the moments of
coexistence equilibrium depends on the transition bias λ and the the stable phenotype distribution but does not affect their shapes
number of phenotypes N. Notably, the stable phenotype qualitatively (Supplementary Fig. 1). We will thus fix N = 3
distribution does not depend on the growth rate as the transition phenotypes in the remainder of the study for illustration.
speed is assumed to be phenotype-independent (see Section
Model for details). Thus, only the transitions and not the growth Transition speed and rate of approach to equilibrium
dynamics decide the phenotype abundances at equilibrium. While the transition bias λ determines the stable phenotype
Additionally, this growth independence shows that the particular distribution, the transition speed c affects how fast the stable
choice of the decline of growth rates from epithelial to phenotype distribution is approached. The model consists of two
mesenchymal phenotypes does not affect the equilibrium processes; logistic growth and transitions, where the parameter c
phenotype distribution. scales the transition propensity relative to the growth rate. When

Published in partnership with the Systems Biology Institute npj Systems Biology and Applications (2023) 48
 S. Shah et al.
4

Fig. 3 Transition speed determines the rate of approach to the stable phenotype distribution. The panels show the approach to the stable
phenotype distribution from the same initial condition, ðx 1 ; x 2 ; x 3 Þ ¼ ð10
1
; 0; 0Þ, for different combinations of transition speed c and transition
bias λ. The transition speed c sets the pace of the transition dynamics relative to the growth dynamics. For any transition bias λ, the time to
reach the stable phenotype distribution decreases as the transition speed increases.

transitions are faster than growth, c > 1, the population first towards the mesenchymal phenotype (Fig. 5c). Growth-
reaches the stable phenotype distribution and then approaches independent treatment targets all phenotypes equally, not
the carrying capacity. When transitions are slower than growth, changing the phenotype distribution directly. The room for
c < 1, the population first reaches the carrying capacity and then growth is filled by epithelial-like cells as they grow fast during
equilibrates to the stable phenotype distribution (Fig. 3). treatment as the tumor is not at carrying capacity (Fig. 5b). Thus,
during growth-independent treatment, the phenotype average
Initial phenotype distribution determines phenotype shifts shifts towards the epithelial phenotype, especially, if transitions
The stability of the coexistence equilibrium induces phenotype are slow compared to growth dynamics (c < 1, Fig. 5d). After
shifts similar to epithelial-mesenchymal and mesenchymal- treatment, the phenotype composition returns to the stable
epithelial transitions. The primary site of a tumor at an early stage distribution that is purely determined by the transition bias (Eq.
is composed mainly of proliferative epithelial cells. Thus, at the (1)). The restoration of the stable phenotype distribution is
primary site, the stable phenotype distribution is approached by a driven by regrowth of the suppressed phenotypes and therefore
relative increase in the abundance of mesenchymal-like cells occurs considerably slower after growth-independent treat-
(Fig. 4, first column). The secondary tumor site is initially founded ment, i.e., in the direction of a frequency increase of
by mesenchymal-like cells due to the selection pressure for mesenchymal cells (Fig. 5d). Such relapses at different rates
invasive properties during dissemination. Therefore, conversely to after distinct treatment types align with our previous findings24
the primary site, the stable equilibrium is approached by a relative and clinical observations (Section 7 in Chitadze et al.40).
increase of epithelial-like cells at the secondary site (Fig. 4, second
column) during metastasis formation. Tumor with high transition speed is vulnerable to treatment
After a perturbation from the coexistence equilibrium, e.g., by Exploring a broader range of values for transition bias and speed,
treatment, the stable phenotype distribution is restored, compen- we find that the effect of growth-dependent and growth-
sating the effect of the perturbation (Fig. 4, columns three to five). independent treatment indeed depends strongly on transition
Hence, the phenotypically plastic tumor cell population can bias and speed (Fig. 6). When a treatment type exerts a particularly
maintain heterogeneity. Our model thus exhibits site- and high mortality on the most abundant tumor phenotype we refer
treatment-specific phenotype shifts that result only from the to the treatment as a phenotype-matched treatment. For instance,
selection during dissemination or treatment. These shifts do not growth-dependent treatment exerts the highest mortality on
require specific microenvironmental factors, although the micro- epithelial tumors, λ → − 1, whereas growth-independent treat-
environment is a crucial component contributing to the initial ment is more phenotype-matched to mesenchymal tumors, λ → 1.
conditions that determine the direction of the shift. We find that one-block treatment schemes are most effective, i.e.
achieve the strongest reduction in tumor burden, when they are
Treatment affects the phenotype distribution of the tumor phenotype-matched (Fig. 6a, b). Applying them for the whole
Next, we tested how the phenotype distribution changes during treatment duration thus results in the highest reduction of tumor
and after growth-dependent and growth-independent treat- burden. Higher transition speeds c cause a faster replenishment of
ment (Fig. 5). Growth-dependent treatment targets proliferating treatment-sensitive phenotypes from less sensitive phenotypes,
cells, thus epithelial cells experience higher treatment-induced entailing a cost of phenotypic plasticity. Thus, treatment efficacy
mortality (Fig. 5a), and the mean of the distribution shifts improves for higher transition speeds for single-block treatments.

npj Systems Biology and Applications (2023) 48 Published in partnership with the Systems Biology Institute
 S. Shah et al.
5

Fig. 4 Plasticity-driven approach to stable phenotype distribution from different initial conditions. Depending on the initial condition, the
approach to the stable phenotype distribution proceeds along different paths (top four rows, the vertical axis represents time). The mean
population phenotype over time is shown in the last row with its asymptote (dashed line). The first column shows the phenotype dynamics of
a plastic tumor at the primary site, where it originates only from epithelial cells. The second column represents the growth of a metastasis
after mainly mesenchymal cells have arrived at the secondary site. The last three columns show the dynamics after a hypothetical intervention
that removes all epithelial, all hybrid, or all mesenchymal cells. In all cases, the tumor approaches the stable phenotype distribution; however,
different initial conditions lead to shifting the average phenotype in different directions.

When the phenotype distribution is known, the adaptive patterns that consist of multiple treatment blocks. In the multi-
treatment scheme can be applied (see Section Treatment types block schemes (see Supplementary Fig. 2 for tumor burden and
and schemes for a detailed description). This treatment scheme mean phenotype dynamics of the two-block scheme), the
selects the treatment type that exerts the higher mortality on the application duration of the optimal treatment is reduced, thus
tumor, given the phenotype distribution at the time of treatment resulting in a less effective treatment of tumors with a high
choice. We assume that initially, the phenotype distribution has transition bias (Fig. 6d, e). Frequent treatment type alterations
reached its equilibrium given by Eq. (1), which is determined result in an intermediate tumor burden reduction across trait
mainly by the transition bias λ. For the one-block adaptive space. Although frequent treatment type alteration is not the
treatment scheme only at the beginning of the treatment phase a most effective strategy, it may be desirable when the character-
treatment type has to be determined. Thus, the treatment type istics of a tumor, such as transition bias λ and transition speed c,
decision in the adaptive one-block treatment scheme depends are unknown (Fig. 6g, h).
only on the transition bias λ. The decision boundary between For the two-block adaptive treatment, the tumor reduction
applying the growth-dependent and the growth-independent pattern is more intricate (Fig. 6f). At λ > ~λ, the mortality exerted by
treatment in terms of transition bias can be found by comparing growth-independent treatment is higher; thus, growth-
the mortality exerted by the two treatmentPtypes. The total independent treatment is chosen as the first treatment type. At
mortality of growth-dependent P treatment N  large transition speed c and high transition bias λ, growth-
i¼1 mD r i x i and
growth-independent treatment N
m x 
are equal when independent treatment is also chosen for the second treatment
PN  
i¼1 I i
block. Only for a transition bias slightly exceeding ~λ and small c, a
i¼1 mD r i x i  mI x i ¼ 0. This is a polynomial in the transition
bias λ with the relevant root ~λ 2 ½1; 1 setting the decision treatment switch occurs, and the tumor burden reduction is
boundary between the two treatment types. For our parameters higher than in the single-block scheme. For λ < ~λ, the first
and assumptions on the growth rates of the phenotypes (Table 1), treatment type is growth-dependent treatment. Only for small c
the switch occurs at ~λ  0:089. Consequently, in the one-block and close to ~λ the treatment type is switched after half the
adaptive treatment (Fig. 6c), the reductions of more epithelial treatment period, which interestingly worsens the treatment
tumors with λ < ~λ are identical to reductions obtained for the outcome compared to a single block of only growth-dependent
treatment. This choice is possible as the decision criterion for
growth-dependent treatment type. More mesenchymal tumors
treatment choice is the instantaneous mortality and does not
with λ > ~λ are affected by one-block adaptive treatment identically
account for the potential future tumor size or the future
to the growth-independent treatment.
phenotype distribution, or mortality integrated over time. Since
the growth-independent treatment exerts no selection pressure
The effect of sequential multi-block treatment schemes against the proliferative phenotype, switching to this treatment
To explore the potential of leveraging the treatment-induced type comes at the risk that soon after the switch, the phenotype
phenotype changes, we also investigated sequential treatment distribution will have shifted to a larger epithelial proportion,

Published in partnership with the Systems Biology Institute npj Systems Biology and Applications (2023) 48
 S. Shah et al.
6

Fig. 5 Growth-dependent and growth-independent treatment can transiently alter the phenotype distribution of a tumor. During
treatment, the tumor burden shrinks (a, b), and the mean phenotype of the tumor changes if the transition speed c is small but restores to the
untreated mean phenotype value after treatment (c, d). Treatment is applied between t = 0 and t = 10. Afterward, regrowth is tracked until
t = 1000. The violet and green horizontal bars indicate growth-dependent and growth-independent treatments. Growth-dependent
treatment shifts the mean of the phenotype distribution towards the mesenchymal phenotype as it exerts higher mortality on epithelial cells.
Conversely, growth-independent treatment shifts the mean towards the epithelial phenotype as epithelial cells compensate for the mortality
by faster growth.

where the growth-independent treatment is less effective than (arrow in Fig. 6b). For multi-block treatment schemes, we find that
the growth-dependent treatment. Increasing the number of increasing transition bias λ generally also increases the reduction
treatment blocks mediates this risk and achieves the highest in tumor burden, albeit this effect depends on the exact pair of
tumor burden reduction for small c and intermediate λ in the transition bias λ and speed c for the adaptive treatment. Adjuvant
adaptive treatment scheme. Different modeling assumptions therapies that increase the transition speed c increase the
(unequal competitiveness of phenotypes or phenotype- treatment effect for a single-block treatment, but for two
dependent transition speed) change the effect of treatment only treatment blocks the outcome is difficult to analyze due to the
quantitatively (see Supplementary Sections Phenotype-dependent dependence on the treatment sequence, transition bias λ and
transition speed, and Unequal competition, Supplementary Figs. 5, speed c.
8).
DISCUSSION
The potential of adjuvant treatment
Cancer cell plasticity increases tumor heterogeneity, facilitates
Modifying the process of phenotype transition with adjuvant metastasis formation, and often hinders treatment approaches41.
therapies promises to alter the phenotype distribution of a tumor, We formulated a mathematical model of phenotypic plasticity that
thus, rendering the tumor more vulnerable to chemo- or motivates a stable heterogeneous phenotype distribution, mainly
immunotherapy. Phenotype transitions may be modulated, determined by the relative propensities of phenotype transitions.
equivalent to changing transition bias λ, by an adjuvant drug We found that the heterogeneity generated by this inherent
that affects the gene regulatory network of epithelial- plasticity can explain the transitions between epithelial and
mesenchymal transitions25,26,29. Figure 6 indicates how modifying mesenchymal phenotypes, which are known to facilitate cancer
the transition bias λ with adjuvant drugs can affect the reduction progression and metastasis4,6,13,14,38,42. Our model provides an
in tumor burden. We find that the tumor burden is reduced more alternative hypothesis for the factors driving such transitions, not
strongly when a single-block growth-dependent treatment is necessarily relying on the microenvironment, yet sensitive to
accompanied by an adjuvant drug that makes the tumor more microenvironmental changes. The stability of a heterogeneous
epithelial by decreasing the transition bias λ (arrow in Fig. 6a). phenotype distribution implies that perturbations to this distribu-
Contrastingly, the effect of the single-block growth-independent tion will be reverted. We discuss the consequences of this
treatment can be enhanced by an adjuvant drug that makes the equilibrium for primary and secondary tumors and their treatment
tumor more mesenchymal by increasing the transition bias λ in this section.

npj Systems Biology and Applications (2023) 48 Published in partnership with the Systems Biology Institute
 S. Shah et al.
7

Fig. 6 Reduction in tumor burden for different sequential treatment schemes. The reduction in tumor burden relative to the carrying
K is indicated by the brightness gradient for different combinations of transition bias λ and transition speed c. Here, X is the sum of
capacity KX
phenotype abundances at the end of treatment duration (N = 3). Darker colors represent a higher reduction and, thus, a better outcome. We
evaluate the effect of splitting the treatment period into multiple treatment blocks (rows) and investigate different treatment schemes with
either predefined or adaptive treatment sequences (columns). The white dashed line indicates a decision boundary ~λ for the adaptive
treatment, which is obtained by comparing the mortality of the two treatment types (see text). Adjuvant therapy can alter the phenotype
transitions, which in our model translates to changes, for example, to the transition bias λ, indicated by the arrows in panels a and b.

Phenotypic plasticity drives phenotype shifts at the primary exacerbate phenotypic changes during the dissemination from
and secondary tumor sites primary to secondary site8. Indeed, we can capture much more
During metastasis formation, one prominent phenotype transition pronounced phenotype changes by assuming different transition
in epithelial cancer cells is the epithelial-mesenchymal transi- biases λ1 ≠ λ2 at primary and secondary sites, resulting in differing
tion14,23. Tumors of epithelial origin gain heterogeneity by stable phenotype distributions.
transitioning to a more mesenchymal phenotype. Selection for
motility and invasiveness during the dissemination of cancer cells Moments of the phenotype distribution guide the choice of
increases the frequency of the mesenchymal phenotype in the phenotype-matched treatment
circulating cancer cells. Therefore, the arriving population at the We investigated growth-dependent and growth-independent
secondary site is mainly constituted of mesenchymal phenotype treatment types, capturing population-level mechanisms of
and gains heterogeneity by recovering more epithelial pheno- chemo- and immunotherapy. We found that the transition bias,
types. Thus, we find that phenotypic plasticity can be responsible and thus the shape of the stable phenotype distribution,
for both the epithelial-mesenchymal transition at the primary site determines the most effective treatment type and sequence of
and the mesenchymal-epithelial transition at the secondary site. different treatment types. Our model exhibits a higher efficiency
However, microenvironmental differences undoubtedly of growth-dependent treatment for a transition bias towards

Published in partnership with the Systems Biology Institute npj Systems Biology and Applications (2023) 48
 S. Shah et al.
8
cancer treatment as the transition bias and, therefore, the stable
Table 1. Symbols, variables and reference parameters.
phenotype distribution can be modulated by an adjuvant therapy.
Symbol Description Value In the context of epithelial-mesenchymal plasticity, the driving
factors of the epithelial-mesenchymal transition and the
xi Abundance of ith tumor phenotype mesenchymal-epithelial transitions are known17,32. Altering the
P
N
X Sum of all phenotype abundances xi concentrations of these driving factors in the tumor microenvir-
i¼1 onment or perturbing the relevant gene regulatory network
N Number of tumor phenotypes 3 would thus translate to changing transition bias λ and speed c
ri Growth rate of ith phenotype r 1  ðr 1  r N Þ i1
N1 which are decisive for the phenotype transition process. For
r1 Growth rate of epithelial cells 1 instance, frequent phenotype switching is costly in rarely
1
changing environments46. Similarly, in our simulations high
rN Growth rate of mesenchymal cells transition speed c improves treatment efficacy. Experimental
5
testing could determine whether indeed faster transition speed
K Carrying capacity of a site 1 poses a cost and can be exploited in clinical settings. Modulation
TEM Transition rate of i → i + 1 c (1 + λ) r1 of transition bias λ has also been shown with TGF-β block-
TME Transition rate of i → i − 1 c (1 − λ) r1 age25,26,28. Therefore, it may be possible to regulate tumor
c Transition speed 1 phenotypes, reduce metastasis, and achieve a better treatment
λ Transition bias 0 response by controlling the transition bias with an adjuvant
treatment. Of note, this reasoning applies also to tumor sites that
mD Growth-dependent treatment intensity 1
have reached carrying capacity. Although the total abundances of
mI Growth-independent treatment intensity ≈ 0.64 such tumors may be static, the phenotype distribution can still
Deviations from these values are reported where applicable. change47. Also in our model, at carrying capacity, the population
growth stops, but the transition dynamics can remain in action,
capturing the potential turnover in tumor composition.
epithelial phenotypes. Growth-independent treatment becomes
more effective for a transition bias towards mesenchymal Driving factors of phenotype transition may serve as a proxy
phenotypes. This connection between treatment effect and the for phenotype abundances
shape of the phenotype distribution can be understood better by
considering the moments of the phenotype distribution. Narrow When measuring the concentrations of driving factors is more
phenotype distributions are well characterized by their mean, feasible than measuring the phenotype abundances, the concen-
which can be used to match treatment type and phenotype. For trations could also be used to characterize tumor heterogeneity
example, chemotherapy will be the most effective when there are building on the decisive effect of such driving factors of the
only proliferative cells. Such a phenotype-matched treatment will phenotype transitions. Hypothetically, in the case of TGF-β driven
perform well for homogeneous tumors in the absence of epithelial-mesenchymal transition, the tumor could be character-
mismatching terminal phenotypes. To implement this strategy in ized by the concentration of TGF-β instead of quantifying the
clinical practice, a thorough characterization of the tumor tissue number of epithelial, hybrid, and mesenchymal cells. More
would be necessary. With the identification of the dominating generally, the state of the microenvironment or the relevant gene
phenotypes a patient may be stratified for phenotype-matched regulatory network can determine the phenotype distribution,
treatment. In-vitro analysis of patient-derived tumor tissues may and thus, be used to inform treatment decisions29.
be used to test and determine the most effective treatment for the Additionally, if available, the third central moment of the
patient. phenotype distribution, i.e., the abundance of mismatched
However, tumors are often heterogeneous, weakening their phenotypes, can characterize the relapse. For example, consider
treatment response43,44. Interestingly, the third central moment of a tumor of epithelial cells with a high third central moment. If this
the phenotype distribution captures the mismatch between a tumor is treated with growth-dependent treatment, then the
treatment type and the abundance of phenotypes that are not relapse will initially be driven by the growth of mesenchymal cells,
matched by the treatment, as they are far off the mean of the eventually, the epithelial cells will take over the population,
phenotype distribution. When matching the treatment to the leading to bi-phasic relapse dynamics24,40. The present model
mean of the phenotype distribution, an increasing absolute third provides a link between the quantification of driving factors of
central moment thus signals the onset of relapse from mis- phenotype transition, equivalent to the transition bias λ, and the
matched phenotypes. In our model, the phenotypic heterogeneity third central moment of the phenotype distribution.
of the stable phenotype distribution is maximal when there is no In this study, we formulated a mathematical model that
transition bias, i.e., λ = 0. The mismatch between targeted quantitatively describes the transition between adjacent pheno-
phenotype and unmatched phenotypes, captured by the third types within a population. Phenotype transitions appear in most
central moment, is largest at intermediate transition bias. forms of unicellular organisms. Our findings, suggesting targeting
Accordingly, our model suggests the existence of a sweet spot phenotype transitions to improve treatment are, hence, applicable
for controlling both heterogeneity and mismatched terminal beyond cancer, i.e., for pathogen-borne diseases such as bacterial
phenotypes at extreme transition biases λ → ± 1, where variance infections. In the current study, we make assumptions that allow
and absolute third central moment are simultaneously minimized. exposition of our findings, e.g. the phenotype transitions are not
adaptive. Also, we do not model the spread between the primary
Driving factors of the phenotype transition control the and secondary cancer sites explicitly. Exchanging some of our
phenotype distribution simple assumptions for realism, we also discuss the case of
Recently, several studies suggested treatment strategies for unequal competition and growth-dependent transition speed in
phenotypically heterogeneous tumors24,35,36,45. These strategies the Supplementary Section Model extensions. The case of unequal
aim to remove or reduce the tumor burden and delay tumor competition can be linked to frequency-dependent interactions
relapse. However, these studies do not target or capitalize on the between tumor phenotypes, which can be captured in the present
source of the heterogeneity. The interaction between treatment model by unequal competition coefficients. There are compelling
and phenotype dynamics can potentially be exploited during studies investigating frequency-dependent interactions48–50 with

npj Systems Biology and Applications (2023) 48 Published in partnership with the Systems Biology Institute
 S. Shah et al.
9
some quantifying them in patient-derived tumors51,52, which could from the growth rates indicates that the particular choice of the
affect treatment outcome and thus, the prognosis of the disease. growth rate decline does not affect our results qualitatively.
To summarize, the present model provides a microenvironment- We assume that the transition rates are equal for all phenotypes
independent explanation for observed phenotype changes during and scaled to the epithelial growth rate r1 with a constant c,
epithelial-mesenchymal and mesenchymal-epithelial transitions. We 0 ≤c < ∞. Thus, transitions faster than proliferation, c > 1, can be
highlight the phenotype transitions in a heterogeneous population interpreted as within-generation phenotypic plasticity, e.g. altered
as a viable clinical target. By targeting the driving factors or the patterns of post-transcriptional regulation, and transitions slower
relevant gene regulatory network of the phenotype transition, we than proliferation c < 1 can be attributed to transgenerational
propose leveraging phenotypic plasticity. Forcing the tumor into a plasticity, such as epigenetic changes. We explore the effect of
less heterogeneous state, this strategy may improve the treatment growth rate-dependent transition rates in the Supplementary
efficacy of chemo- and immunotherapies. Supporting the estab- Section Phenotype-dependent transition speed.
lished idea of personalized treatments, our approach thus discusses To account for a potential bias in the direction of plastic
a practical option to bring the concept into clinical practice. transitions, we introduce the transition bias λ, − 1 ≤ λ ≤ 1, which
determines the probability of a cell switching to the adjacent
epithelial or mesenchymal phenotype. A positive transition bias
METHODS λ > 0 implies a bias towards becoming more mesenchymal, λ < 0
We investigate the onset of metastasis formation mediated by the results in a higher probability of switching to a more epithelial
plasticity between epithelial and mesenchymal phenotypes. We phenotype. With these assumptions, the transition rates become
use an ecology-inspired approach to identify tumor phenotypes TEM = c (1 + λ) r1 and TME = c (1 − λ) r1. The parameters of our
with species in an ecosystem, drawing on the usefulness of model are given in Table 1.
ecological concepts for understanding cancer biology53,54. We use
the competitive Lotka–Volterra system55,56 with additional terms Treatment types and schemes
that capture transitions between phenotypes. We go beyond To investigate the effect of drug treatment on phenotypically
earlier approaches that featured two cell types35,36 and account plastic tumor cell populations, we apply two different treatment
for the recently established prominence of hybrid phenotypes types. The mortality exerted on the cells of a particular phenotype
characterized by both partial epithelial and mesenchymal by the two treatment types either does or does not depend on the
markers15. This results in a multi-compartment model similar to phenotype’s growth rate. The death rate exerted by growth-
those presented in Zhou et al.33 and Raatz et al.24, but with added dependent treatment is mD ri, and the death rate by growth-
competition and more flexibility in specifying the transition independent treatment is mI. To better compare the phenotype-
process. Using this model, we examine the epithelial- distribution-dependent treatment effects, we chose mI such that
mesenchymal plasticity and characterize how different properties the reduction of tumor cell abundance at the end of the treatment
of the plastic transition shape the resulting tumor heterogeneity. duration is identical for both treatment types. To find mI, we use
c = 1, λ = 0, mD = 1, and numerically minimize the difference
between the total abundances X at the end of treatment duration.
Model
Throughout this study, treatment is always applied for a fixed
We represent the phenotypic heterogeneity as a collection of N duration of 10 time units.
different tumor phenotypes at different sites of cancer (Eq. (2)). These assumptions result in a system of ordinary differential
These phenotypes range from the most epithelial phenotype, E, equations describing the change in tumor phenotype abundance,
over N − 2 hybrid phenotypes H, to the most mesenchymal
phenotype, M. We assume that all tumor phenotypes share  
dt ¼ r1 x1 1  XK þ c r1 ½ð1  λÞ x2  ð1 þ λÞ x1   ðmD r1 þ mI Þ x1
dx1
resources at a given cancer site and that the tumor population
follows logistic growth57. Additionally, we assume that the tumor ..
.  
cells are phenotypically plastic, i.e., they can change to an X
adjacent, more epithelial-like, or mesenchymal-like phenotype in
dxi
¼ ri x i 1  þ c r1 ½ð1 þ λÞ xi1 þ ð1  λÞ xiþ1  2 xi   ðmD ri þ mI Þ xi :
dt K |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}
bidirectional phenotype transitions. |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} transitions treatments
logistic growth
We track the phenotype abundance xi over time for each tumor ..
phenotype i, i = 1, 2, …, N. Each cell of tumor phenotype i grows .
 
intrinsically at rate ri. The cells switch to the next functionally more dxN
dt ¼ rN xN 1  XK þ c r1 ½ð1 þ λÞ xN1  ð1  λÞ xN   ðmD rN þ mI Þ xN
mesenchymal phenotype, i → i + 1, at rate TEM, and to the next
(2)
functionally more epithelial type, i → i − 1, at rate TME. The positive
transition terms represent additions to a phenotype from adjacent Different drug administration strategies are possible and of clinical
phenotypes. The negative transition terms account for the losses interest. However, in current practice, mostly monotherapies or
from transitions out of the current phenotype. Resource scarcity in combined therapies are used, often leading to resistance due to
the microenvironment induces a density-dependent, logistic static selection pressures. This motivates us to investigate
P sequential treatment schemes. The fixed treatment scheme either
competition captured by r i XK , where X ¼ Ni¼1 x i is the total
population abundance and K is the carrying capacity. We assume treats with the same drug throughout the treatment duration or
alternates between blocks of different treatment types. The single-
that all phenotypes are equally competitive for resources, i.e., the
block treatment scheme applies either growth-dependent or
competition term is equal for all phenotypes. Relaxing this
growth-independent treatment. In the multi-block treatment
assumption does not affect our results qualitatively and is
schemes, we alternate between treatment types and investigate
explored in the Supplementary Section Unequal competition.
both cases of first treating with growth-dependent or first treating
Further, we assume linearly decreasing growth rates from
with growth-independent treatment. For the adaptive treatment
epithelial to mesenchymal phenotypes, reflecting the higher scheme, at the start of each treatment block, the treatment type is
proliferation rates of epithelial cells. We scale all growth rates chosen that achieves the higher instantaneous reduction of the
relative to the epithelial growth rate r1, and express the growth whole tumor cell population for the given phenotype distribu-
rate for the ith phenotype as r i ¼ r 1  ðr 1  r N Þ N1
i1
. Note that all tion24. The cancer cell mortality for growth-dependent treatment
the growth rates are equally distributed in the interval [r1, rN]. The P
is Ni¼1 mD r i x i and the mortality exerted by growth-independent
independence of the equilibrium phenotype distribution (Eq. (1))

Published in partnership with the Systems Biology Institute npj Systems Biology and Applications (2023) 48
 S. Shah et al.
10
PN
treatment is 16. Hong, D. et al. Epithelial-to-mesenchymal transition and cancer stem cells con-
i¼1 mI x i . The comparison of these two terms
determines the decision for the treatment type. tribute to breast cancer heterogeneity. J. Cell. Physiol. 233, 9136–9144 (2018).
17. Dongre, A. & Weinberg, R. A. New insights into the mechanisms of
epithelial–mesenchymal transition and implications for cancer. Nat. Rev. Mol. Cell
Analysis and implementation Biol. 20, 69–84 (2019).
We performed linear and Lyapunov stability analysis to find the 18. Jordan, N. V., Johnson, G. L. & Abell, A. N. Tracking the intermediate stages of
equilibria of Eqs. (2) and establish their stability58. The results of epithelial-mesenchymal transition in epithelial stem cells and cancer. Cell Cycle
10, 2865–2873 (2011).
this analysis are presented in the Supplementary Section Stability
19. Huang, R. Y.-J. et al. An EMT spectrum defines an anoikis-resistant and spher-
of the coexistence equilibrium. For N = 2 and N = 3 stability was oidogenic intermediate mesenchymal state that is sensitive to e-cadherin
confirmed computationally (Supplementary Fig. 3), for N > 3 we restoration by a src-kinase inhibitor, saracatinib (AZD0530). Cell Death Dis. 4,
show stability using Lyapunov stability analysis. Temporal e915–e915 (2013).
dynamics of the model (Eqs. (2)) were obtained by numerical 20. Jolly, M. K. et al. Hybrid epithelial/mesenchymal phenotypes promote metastasis and
integration using the solve_ivp function from the Scipy therapy resistance across carcinomas. Pharmacol. Therapeutics 194, 161–184 (2019).
library59 in Python (version 3.960). Numpy61 and Matplotlib62 were 21. Selvaggio, G. et al. Hybrid epithelial–mesenchymal phenotypes are controlled by
used for computation and plotting. microenvironmental factors. Cancer Res. 80, 2407–2420 (2020).
22. Goetz, H., Melendez-Alvarez, J. R., Chen, L. & Tian, X.-J. A plausible accelerating
function of intermediate states in cancer metastasis. PLoS Comput. Biol. 16,
Reporting summary e1007682 (2020).
Further information on research design is available in the Nature 23. Nieto, M. A., Huang, R. Y.-J., Jackson, R. A. & Thiery, J. P. EMT: 2016. Cell 166, 21–45
Research Reporting Summary linked to this article. (2016).
24. Raatz, M., Shah, S., Chitadze, G., Brüggemann, M. & Traulsen, A. The impact of
phenotypic heterogeneity of tumour cells on treatment and relapse dynamics.
DATA AVAILABILITY PLoS Comput. Biol. 17, e1008702 (2021).
25. Liu, J. et al. TGF-β blockade improves the distribution and efficacy of therapeutics
The figures and dataset generated in this study have been deposited in a zenodo
in breast carcinoma by normalizing the tumor stroma. Proc. Natl Acad. Sci. 109,
repository with https://doi.org/10.5281/zenodo.7989753.
16618–16623 (2012).
26. Pickup, M., Novitskiy, S. & Moses, H. L. The roles of TGFβ in the tumour micro-
environment. Nat. Rev. Cancer 13, 788–799 (2013).
CODE AVAILABILITY 27. Wang, W., Poe, D., Yang, Y., Hyatt, T. & Xing, J. Epithelial-to-mesenchymal tran-
The relevant code to generate the figures and dataset can be found in a zenodo sition proceeds through directional destabilization of multidimensional attractor.
repository with https://doi.org/10.5281/zenodo.7989748. eLife 11, e74866 (2022).
28. Lee, H.-M., Lee, H.-J. & Chang, J.-E. Inflammatory cytokine: an attractive target for
Received: 9 February 2023; Accepted: 15 September 2023; cancer treatment. Biomedicines 10, 2116 (2022).
29. Hari, K., Ullanat, V., Balasubramanian, A., Gopalan, A. & Jolly, M. K. Landscape of
epithelial–mesenchymal plasticity as an emergent property of coordinated teams
in regulatory networks. eLife 11, e76535 (2022).
30. Lu, M., Jolly, M. K., Levine, H., Onuchic, J. N. & Ben-Jacob, E. MicroRNA-based
regulation of epithelial–hybrid–mesenchymal fate determination. Proc. Natl Acad.
REFERENCES Sci. 110, 18144–18149 (2013).
1. Dillekås, H., Rogers, M. S. & Straume, O. Are 90% of deaths from cancer caused by 31. Deshmukh, A. P. et al. Identification of EMT signaling cross-talk and gene reg-
metastases? Cancer Med. 8, 5574–5576 (2019). ulatory networks by single-cell RNA sequencing. Proc. Natl Acad. Sci. 118,
2. Siegel, R. L., Miller, K. D. & Jemal, A. Cancer statistics, 2020. Cancer J. Clin. 70, 7–30 e2102050118 (2021).
(2020). 32. Subbalakshmi, A. R. et al. KLF4 induces mesenchymal–epithelial transition (MET)
3. Sung, H. et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence by suppressing multiple EMT-inducing transcription factors. Cancers 13, 5135
and Mortality Worldwide for 36 Cancers in 185 Countries. Cancer J. Clin. 71, (2021).
209–249 (2021). 33. Zhou, D., Luo, Y., Dingli, D. & Traulsen, A. The invasion of de-differentiating cancer
4. Gupta, G. P. & Massagué, J. Cancer metastasis: building a framework. Cell 127, cells into hierarchical tissues. PLoS Comput. Biol. 15, e1007167 (2019).
679–695 (2006). 34. Gupta, P. B. et al. Stochastic state transitions give rise to phenotypic equilibrium
5. Weinberg, R. A. The Biology of Cancer 2nd edn (Garland Science, Taylor & Francis in populations of cancer cells. Cell 146, 633–644 (2011).
Group, 2014). ISBN 978-0-8153-4219-9 978-0-8153-4220-5. 35. Li, X. & Thirumalai, D. A mathematical model for phenotypic heterogeneity in
6. Turajlic, S. & Swanton, C. Metastasis as an evolutionary process. Science 352, breast cancer with implications for therapeutic strategies. J. Roy. Soc. Interface 19,
169–175 (2016). 20210803 (2022).
7. Shlyakhtina, Y., Moran, K. L. & Portal, M. M. Genetic and non-genetic mechanisms 36. Cassidy, T., Nichol, D., Robertson-Tessi, M., Craig, M. & Anderson, A. R. A. The role
underlying cancer evolution. Cancers 13, 1380 (2021). of memory in non-genetic inheritance and its impact on cancer treatment
8. Group Young Researchers in Inflammatory Carcinogenesis, Wandmacher, A. M., resistance. PLoS Comput. Biol. 17, e1009348 (2021).
Mehdorn, A.-S. & Sebens, S. The heterogeneity of the tumor microenvironment as 37. Ma, B., Wells, A. & Clark, A. M. The pan-therapeutic resistance of disseminated
essential determinant of development, progression and therapy response of tumor cells: Role of phenotypic plasticity and the metastatic microenvironment.
pancreatic cancer. Cancers 13, 4932 (2021). Semin. Cancer Biol. 60, 138–147 (2020).
9. Caiado, F., Silva-Santos, B. & Norell, H. Intra-tumour heterogeneity—going 38. Esposito, M., Ganesan, S. & Kang, Y. Emerging strategies for treating metastasis.
beyond genetics. FEBS J. 283, 2245–2258 (2016). Nat. Cancer 2, 258–270 (2021).
10. West-Eberhard, M. J. Phenotypic plasticity and the origins of diversity. Annu. Rev. 39. Franssen, L. C. & Chaplain, M. A. J. A mathematical multi-organ model for bidir-
Ecol. Systematics 20, 249–278 (1989). ectional epithelial–mesenchymal transitions in the metastatic spread of cancer.
11. Fischer, B. B. et al. Phenotypic plasticity influences the eco-evolutionary dynamics IMA J. Appl. Math. 85, 724–761 (2020).
of a predator–prey system. Ecology 95, 3080–3092 (2014). 40. Chitadze, G., Laqua, A., Lettau, M., Baldus, C. D. & Brüggemann, M. Bispecific
12. Giese, A. et al. Dichotomy of astrocytoma migration and proliferation. Int. J. antibodies in acute lymphoblastic leukemia therapy. Expert Rev. Hematol. 13,
Cancer 67, 275–282 (1996). 1211–1233 (2020).
13. Zheng, X. et al. Epithelial-to-mesenchymal transition is dispensable for metas- 41. Deshmukh, S. & Saini, S. Phenotypic heterogeneity in tumor progression, and its
tasis but induces chemoresistance in pancreatic cancer. Nature 527, 525–530 possible role in the onset of cancer. Front. Genet. 11, 604528 (2020).
(2015). 42. Fares, J., Fares, M. Y., Khachfe, H. H., Salhab, H. A. & Fares, Y. Molecular principles
14. Jolly, M. K., Ware, K. E., Gilja, S., Somarelli, J. A. & Levine, H. EMT and MET: of metastasis: a Hallmark of cancer revisited. Signal Transduct. Targeted Ther. 5, 28
Necessary or permissive for metastasis? Mol. Oncol. 11, 755–769 (2017). (2020).
15. Jolly, M. K., Mani, S. A. & Levine, H. Hybrid epithelial/mesenchymal phenotype(s): 43. Greenbaum, A. et al. Tumor heterogeneity as a predictor of response to
The ‘fittest’ for metastasis? Biochim. Biophys. Acta (BBA)—Rev. Cancer 1870, neoadjuvant chemotherapy in locally advanced rectal cancer. Clin. Colorectal
151–157 (2018). Cancer 18, 102–109 (2019).

npj Systems Biology and Applications (2023) 48 Published in partnership with the Systems Biology Institute
 S. Shah et al.
11
44. McDonald, K.-A. et al. Tumor heterogeneity correlates with less immune response the Deutsche Krebshilfe (AZ 70112935) to Su.Se. Sa.Sh. is grateful to Dr. Florence
and worse survival in breast cancer patients. Ann. Surg. Oncol. 26, 2191–2199 Bansept for suggesting alternate parametrization and to Saptarshi Pal for the
(2019). discussions on stability.
45. Viossat, Y. & Noble, R. A theoretical analysis of tumour containment. Nat. Ecol.
Evol. 5, 826–835 (2021).
46. Acar, M., Mettetal, J. T. & van Oudenaarden, A. Stochastic switching as a survival AUTHOR CONTRIBUTIONS
strategy in fluctuating environments. Nat. Genet. 40, 471–475 (2008). Conceptualization: Sa.Sh., M.R., A.T., Su.Se., L.M.P., Methodology: Sa.Sh., M.R., A.T.,
47. Raatz, M. & Traulsen, A. Promoting extinction or minimizing growth? The impact Software: Sa.Sh., Investigation: Sa.Sh., M.R., S.G., A.T., Visualization: Sa.Sh., Writing—
of treatment on trait trajectories in evolving populations. Evolution qpad042, original draft: Sa.Sh., M.R., Writing—review and editing: Sa.Sh., M.R., A.T., S.G. Su.Se.,
March (2023). L.M.P., Supervision and project administration: M.R., A.T.
48. Tomlinson, I. P. M. Game-theory models of interactions between tumour cells.
Eur. J. Cancer 33, 1495–1500 (1997).
49. Tomlinson, I. & Bodmer, W. F. Modelling the consequences of interactions
FUNDING
between tumour cells. Br. J. Cancer 75, 157–160 (1997).
50. Böttcher, M. A. et al. Modeling treatment-dependent glioma growth including a Open Access funding enabled and organized by Projekt DEAL.
dormant tumor cell subpopulation. BMC Cancer 18, 376 (2018).
51. Freischel, A. R. et al. Frequency-dependent interactions determine outcome of
competition between two breast cancer cell lines. Sci. Rep. 11, 4908 (2021). COMPETING INTERESTS
52. Kaznatcheev, A., Peacock, J., Basanta, D., Marusyk, A. & Scott, J. G. Fibroblasts and The authors declare no competing interests.
alectinib switch the evolutionary games played by non-small cell lung cancer.
Nat. Ecol. Evol. 3, 450–456 (2019).
53. Basanta, D. & Anderson, A. R. A. Exploiting ecological principles to better ADDITIONAL INFORMATION
understand cancer progression and treatment. Interface Focus 3, 20130020 Supplementary information The online version contains supplementary material
(2013). available at https://doi.org/10.1038/s41540-023-00309-1.
54. Korolev, K. S., Xavier, J. B. & Gore, J. Turning ecology and evolution against cancer.
Nat. Rev. Cancer 14, 371–380 (2014). Correspondence and requests for materials should be addressed to Saumil Shah.
55. Smale, S. On the differential equations of species in competition. J. Math. Biol. 3,
5–7 (1976). Reprints and permission information is available at http://www.nature.com/
56. Hofbauer, J. & Sigmund, K. Evolutionary Games and Population Dynamics 1st edn reprints
(Cambridge University Press, 1998) ISBN 978-0-521-62365-0 978-0-521-62570-8
978-1-139-17317-9. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims
57. Vaidya, V. G. & Alexandro, F. J. Evaluation of some mathematical models for in published maps and institutional affiliations.
tumor growth. Int. J. Bio-Med. Comput. 13, 19–35 (1982).
58. Strogatz, S. H. Nonlinear Dynamics and Chaos 0th edn (CRC Press, 2018). ISBN 978-
0-429-96111-3.
59. Virtanen, P. et al. SciPy 1.0: Fundamental algorithms for scientific computing in Open Access This article is licensed under a Creative Commons
Python. Nat. Methods 17, 261–272 (2020). Attribution 4.0 International License, which permits use, sharing,
60. van Rossum, G. and Drake, F. L.The Python Language Reference. Number Pt. 2 in adaptation, distribution and reproduction in any medium or format, as long as you give
Python Documentation Manual / Guido van Rossum; Fred L. Drake [Ed.]. Python appropriate credit to the original author(s) and the source, provide a link to the Creative
Software Foundation, Hampton, NH, release 3.0.1 [repr.] edition, (2010). ISBN 978- Commons license, and indicate if changes were made. The images or other third party
1-4414-1269-0. material in this article are included in the article’s Creative Commons license, unless
61. Harris, C. R. et al. Array programming with NumPy. Nature 585, 357–362 (2020). indicated otherwise in a credit line to the material. If material is not included in the
62. Hunter, J. D. Matplotlib: A 2D graphics environment. Comput. Sci. Eng. 9, 90–95 article’s Creative Commons license and your intended use is not permitted by statutory
(2007). regulation or exceeds the permitted use, you will need to obtain permission directly
from the copyright holder. To view a copy of this license, visit http://
creativecommons.org/licenses/by/4.0/.
ACKNOWLEDGEMENTS
All authors acknowledge funding by Deutsche Forschungsgemeinschaft through the
Research Training Group “Translational Evolutionary Research” (TransEvo) (Project © The Author(s) 2023
number 400993799, https://gepris.dfg.de/gepris/projekt/400993799) and funding by

Published in partnership with the Systems Biology Institute npj Systems Biology and Applications (2023) 48