Stochastic modeling of phenotypic switching and chemoresistance in cancer cell populations

TL;DR

This paper develops stochastic models to understand how phenotypic switching, especially related to EMT, contributes to drug resistance in cancer cells, providing analytical tools to analyze heterogeneity and tumor response.

Contribution

It introduces a coarse-grained stochastic model linking phenotypic switching and drug resistance, supported by experimental observations, and derives analytical results for cell heterogeneity.

Findings

Analytical probability distributions for phenotypic states over time

Insights into switching rates and heterogeneity in cancer cell populations

Framework for understanding tumor adaptation to chemotherapy

Abstract

Phenotypic heterogeneity in cancer cells is widely observed and is often linked to drug resistance. In several cases, such heterogeneity in drug sensitivity of tumors is driven by stochastic and reversible acquisition of a drug tolerant phenotype by individual cells even in an isogenic population. Accumulating evidence further suggests that cell-fate transitions such as the epithelial to mesenchymal transition (EMT) are associated with drug resistance. In this study, we analyze stochastic models of phenotypic switching to provide a framework for analyzing cell-fate transitions such as EMT as a source of phenotypic variability in drug sensitivity. Motivated by our cell-culture based experimental observations connecting phenotypic switching in EMT and drug resistance, we analyze a coarse-grained model of phenotypic switching between two states in the presence of cytotoxic stress from chemotherapy. We derive analytical results for time-dependent probability distributions that provide insights into the rates of phenotypic switching and characterize initial phenotypic heterogeneity of cancer cells. The results obtained can also shed light on fundamental questions relating to adaptation and selection scenarios in tumor response to cytotoxic therapy.

Figures

Equations

$$\frac{\partial P(E,M,t)}{\partial t}$$

$$g(z_1,z_2,t|p_0)=1-(P_E+P_M)+P_E{z}_1+P_M{z}_2.$$

$$G(z_1,z_2,t)={\int }_{p_0=0}^{p_0=1}d{p}_0\rho (p_0)g(z_1,z_2,t|p_0{)}^{N_0},$$

$$⟨N⟩/N_0$$

$${\gamma }_0=k_{EM}+k_{ME}+{\mu }_E+{\mu }_M,{\alpha }_0=\sqrt{{\gamma }_0^2-4\left(k_{ME}({\mu }_E-{\mu }_M)+({\gamma }_0-{\mu }_M){\mu }_M\right)}.$$

$${\sigma }_N^2={\sigma }_E^2+{\sigma }_M^2+2C_{EM},$$

$$F=1-\frac{⟨N⟩}{N_0}+\frac{N_0(N_0-1)}{⟨N⟩}{\left[\left(\frac{({\mu }_E-{\mu }_M)}{{\alpha }_0}\right)\left(\exp \left(-\frac{t}{2}({\gamma }_0-{\alpha }_0)\right)-\exp \left(-\frac{t}{2}({\gamma }_0+{\alpha }_0)\right)\right)\right]}^2{\sigma }_{p_0}^2.$$

$$\mathcal{F}=\left(\frac{{\sigma }_{p_0}({\mu }_E-{\mu }_M)}{{\alpha }_0}\right)\left(\exp \left(\frac{t}{2}({\gamma }_0-{\alpha }_0)\right)-\exp \left(\frac{t}{2}({\gamma }_0+{\alpha }_0)\right)\right),$$

$$\mathcal{F}=\exp \left[\frac{\ln \left(\left(F-1+\frac{⟨N⟩}{N_0}\right)\frac{⟨N⟩}{N_0(N_0-1)}\right)}{2}\right].$$

$$p_0$$

$${\sigma }_{p_0}^2$$

$$⟨N⟩$$

$$\gamma$$

$$\frac{\partial ⟨E^2⟩}{\partial t}$$

$$\frac{\partial ⟨M^2⟩}{\partial t}$$

$$\frac{\partial ⟨ME⟩}{\partial t}$$

$$\rho (p_0)=\frac{\Gamma (\alpha +\beta )}{\Gamma (\alpha )\Gamma (\beta )}{p}_0^{\alpha -1}(1-p_0{)}^{\beta -1}$$

$$\alpha =\frac{⟨p_0⟩\left(⟨p_0⟩(1-⟨p_0⟩)-{\sigma }_{p_0}^2\right)}{{\sigma }_{p_0}^2},\beta =\frac{(⟨p_0⟩-1)\left[⟨p_0{⟩}^2-⟨p_0⟩+{\sigma }_{p_0}^2\right]}{{\sigma }_{p_0}^2}.$$

$$g(z_1,z_2,t|p_0)=P(0,0,t|p_0)+z_{1}P(1,0,t|p_0)+z_{2}P(0,1,t|p_0).$$

$$\frac{d{P}_E(t)}{dt}$$

$$\frac{d{P}_M(t)}{dt}$$

$$P_E$$

$$P_M$$

$$G(z_1,z_2,t)=\sum _{E=0}^{\infty }\sum _{M=0}^{\infty }{z}_1^E{z}_2^{M}P(E,M,t)={\int }_{p_0=0}^{p_0=1}d{p}_0\rho (p_0){\left[g(z_1,z_2,t)\right]}^{N_0},$$

$$⟨E⟩$$

$$⟨M⟩$$

$$Q_0$$

$$Q_1$$

$$S_0$$

$$S_1$$

$$⟨N⟩=N_0\left[Q_0+S_0+(Q_1+S_1)⟨p_0⟩\right].$$

$${\sigma }_E^2={\frac{d^{2}G}{d{z}_1^2}|}_{1,1}-{\left({\frac{dG}{d{z}_1}|}_{1,1}\right)}^2+{\frac{dG}{d{z}_1}|}_{1,1},$$

$${\sigma }_M^2={\frac{d^{2}G}{d{z}_2^2}|}_{1,1}-{\left({\frac{dG}{d{z}_2}|}_{1,1}\right)}^2+{\frac{dG}{d{z}_2}|}_{1,1},$$

$$F_E$$

$$F_M$$

$$F_E$$

Full text

I Introduction

Acquisition of drug resistance constitutes a major challenge in cancer therapy Housman et al. (2014); Gottesman (2002); Chisholm et al. (2016); Bozic and Nowak (2017); Pogrebniak and Curtis (2018); Nikolaou et al. (2018); Salgia and Kulkarni (2018); Dagogo-Jack and Shaw (2018); Zhou et al. (2009); Zahreddine and Borden (2013); Holohan et al. (2013); Garraway and Jänne (2012); Shaffer et al. (2017). Therapeutic agents (with widely varying biochemical mechanisms) often exhibit a common pattern of providing an initial reduction in tumor burden followed by recurrence of therapeutically resistant disease with more aggressive progression Nikolaou et al. (2018); Holohan et al. (2013); Garraway and Jänne (2012); Chabner and Roberts (2005). Tumor recurrence, which is a major obstacle for cancer cure, is primarily associated with the survival and growth of cell phenotypes that are resistant to chemotherapy Gatenby and Brown (2018); Gallaher et al. (2018); Castorina et al. (2009). Therefore, in order to develop new strategies for the effective treatment of human cancers, a quantitative understanding of the underlying processes leading to drug resistance is essential.

Cellular phenotypic heterogeneity is widely observed in many cancers Pardal et al. (2003); Chisholm et al. (2016); Pogrebniak and Curtis (2018); Dagogo-Jack and Shaw (2018); Meacham and Morrison (2013); Marusyk et al. (2012) as a tumor is often composed of multiple subpopulations Gupta et al. (2009); Zhou et al. (2014) that show different responses to chemotherapy Zhou et al. (2009). In particular, cellular phenotypes that are not sensitive to drugs survive the treatment and can drive drug resistance. As the underlying processes that can lead to the emergence of resistant cells are often stochastic, tumors may locally contain varying numbers of resistant cells. Therefore, quantifying the statistics of drug resistant cells in a tumor is important for effective therapy. Specifically, we are interested in studying population heterogeneity at the start of therapy and aim to address an important issue of therapeutic importance, namely, how to quantify randomness in the numbers of resistant cells prior to drug treatment.

In analyzing population heterogeneity in tumors, a basic question that arises is: How are cell phenotypes that confer survival advantage in the presence of chemotherapeutic drugs generated? A common explanation for the emergence of such phenotypes revolves around Darwinian selection of pre-existing cellular heterogeneity that arises due to random genetic mutations Nowell (1976); Sottoriva et al. (2013); Burrell et al. (2013). However, the fact that resistant cells switch reversibly to sensitive cells, and that resistant cells often appear on short time intervals (hours to few days), starting from clonal populations, suggests that non-genetic factors play a major role in the generation of phenotypic heterogeneity Pisco et al. (2013); Pisco and Huang (2015); Pogrebniak and Curtis (2018); Salgia and Kulkarni (2018); Brown et al. (2014); Shaffer et al. (2017); Su et al. (2017); Inde and Dixon (2018). Such non-genetic phenotypic heterogeneity can arise due to multistability in the underlying gene expression dynamics Chang et al. (2008); Huang et al. (2005) and noise in gene expressionKaern et al. (2005).

These observations suggest that there are two distinct, though not mutually exclusive, mechanisms for the onset of drug resistance in cancer cells: 1) cell phenotypes that are resistant to chemotherapy pre-exist in the tumor prior to treatment and are selected for during the treatment, and 2) cells are induced to develop or acquire resistance due to treatment. In the literature, the former scenario is termed selection while later corresponds to adaptation. This adaptation-selection scenario was first explored in the famous Luria-Delbrück experiments Luria and Delbrück (1943) to understand the mechanism of bacterial resistance to bacteriophage infections. The corresponding analysis gave rise to the celebrated fluctuation test which is also used to estimate mutation rates in bacteria. It is important to note that, while in the Luria-Delbrück case phenotypic changes are driven by genetic mutations and thus an irreversible process, in our study, we are considering phenotypic changes that are reversible. Besides reversible phenotypic switching, it is important to consider intrinsic stochasticity in the underlying processes and to characterize cellular heterogeneity as highlighted by previous studies focusing on modeling drug resistance in cancer Kessler et al. (2014); Komarova (2006).

In consideration of cellular mechanisms likely to be associated with drug resistance, the epithelial-mesenchymal transition (EMT) emerges as a logical candidate. EMT is a conserved cellular program that enables cells of epithelial lineage to transiently acquire traits of mesenchymal cells, including reversible loss of adherens junctions and gain of proteins associated with enhanced motility, adhesion to extracellular substrates and remodeling of the extracellular matrix Kalluri and Weinberg (2009); Lamouille et al. (2014). In cancer cells, this ability to reversibly adopt a more motile phenotype has been linked to tumor invasion and metastasis Heerboth et al. (2015); Yang and Weinberg (2008); Zhang and Weinberg (2018) but, more importantly for this study, EMT is also directly linked to chemotherapy resistance and cancer stem cell (CSC) properties Thiery et al. (2009); Singh and Settleman (2010). In the context of this background, experimental studies described herein focus on established markers of epithelial and mesenchymal phenotype in relation to chemotherapy response, which in this report involves pancreatic ductal adenocarcinoma (PDAC) cells. While recognizing that EMT is more likely a spectrum of intermediate states Lu et al. (2013); Jolly et al. (2015, 2016); Hong et al. (2015); Li and Balazsi (2018), the strong correlation in phenotype and drug response reported here motivates the adoption of two coarse-grained states to be used in the model development. Specifically, a relatively drug-sensitive state with more pronounced epithelial characteristics (E); and a drug-resistant state with increased mesenchymal characteristics (M). In the following sections, we will consider these phenotypes to form the basis of a two-state model of the dynamics of phenotypic switching and associated survival of cancer cells under cytotoxic stressPisco et al. (2013); Zhou et al. (2014).

The paper is organized as follows. In Section II, a set of motivating experiments is described, in which drug resistance is evaluated as a determinant of phenotype in pancreatic cancer cells in vitro, and conversely, phenotype as a determinant of drug response in the same cells. Motivated by these studies, a two phenotype switching model is described in Sec III. An analytic approach to quantify population heterogeneity at the start of therapy is presented in the Sec IV. Then, a protocol for estimating switching parameters is presented in Sec V. In Sec VI, we present an approximate approach for characterizing the probability distribution of the fraction of resistant cells in a population followed by conclusions presented in Sec VII.

II Evaluation of drug resistance and phenotype in cell culture studies

We first sought to compare phenotypic traits in naive and drug-resistant pancreatic ductal adenocarcinoma (PDAC) cells. PANC1 cells (a quasimesnchymal PDAC cell line Collisson et al. (2011)) were exposed to increasing doses of oxaliplatin chemotherapy over successive passages until resistant cells were stable through multiple passages and cryopreservation. As shown in Figure 1 (upper panels), acquisition of chemoresistance leads to a marked change in phenotype from naive cells displaying characteristically epithelial adherens junctions, to drug-resistant cells with highly branched morphology, no evident E-cadherin, and marked increase in cytoskeletal vimentin (IF quantification, upper right). This pattern of changes in E-cadherin and vimentin expression are classic and well-established markers of EMT. We further examined the reciprocal scenario, in which the same parental cells were directly induced Kalluri and Weinberg (2009)to undergo EMT via administration of exogenous TGF-$\beta$ (Figure 1, lower panels). The resultant phenotype is strikingly similar to that of our drug-resistant cells and importantly, exhibits resistance to chemotherapy similar to when resistance was acquired directly through drug exposure. Collectively these results display a symmetry in that acquisition of drug resistance in epithelial cancer cells leads to increase in mesenchymal characteristics, while direct transition from epithelial to mesenchymal phenotype leads to drug resistance.

Methods

Cell culture and reagents: PANC1 cells were obtained from the American Type Culture Collection (Manassas, VA), and grown in T-75 cell culture flasks according to ATCC guidelines. DMEM medium (HyClone; Logan, UT) was supplemented with 10% FBS (HyClone; Logan, UT), 100 IU/mL penicillin and 1% streptomycin (HyClone; Logan, UT), and 0.5 ug/mL Amphotericin B (Corning; Corning, NY). The drug-resistant subline, PANC1-OR was generated as described previously Cramer et al. (2017). Briefly, increasing concentrations of oxaliplatin were added to each cell type in regular media over the course of approximately 25 passages until a stable proliferative phenotype without chemotherapy was observed and maintained following cryopreservation and confirmed by comparative dose response and measurement of a statistically significant increase in IC50.

Immunofluorescence sample preparation and imaging: Formaldehyde-fixed cells in optical-bottom multiwell plates were incubated overnight at $4^{\circ}$C with primary antibodies against e-cadherin and vimentin (Cell Signaling EMT Duplex; Danvers, MA). After washing with PBS, cells were incubated for 1 hour with mouse or rabbit Alexa Fluor secondary antibodies (Cell Signaling; Danvers, MA). Cells were mounted with ProLong Gold Antifade reagent containing DAPI (ThermoFisher Scientific Molecular Probes; Waltham, MA) and imaged after 24 hours using a Zeiss LSM 880 confocal microscope with the same detector settings and excitation laser power settings across groups. Images were analyzed using custom Matlab scripts where fluorescent signal for each protein was normalized to the number of cells based on DAPI-stained nuclei.

Therapeutic response assessment: In sample wells receiving chemotherapy treatment, oxaliplatin (Selleck Chemical; Houston, TX) was added to the media at doses ranging from 0.1 to 500 $\mu$M for 48 hours. In experiments where EMT was induced via TGF-beta, 10 ng/mL human recombinant TGF-beta (Gibco, Thermo Fisher Scientific) in 1% FBS DMEM was added to designated wells for 48 hours and respective comparison groups were also grown in 1% FBS for the same duration. In all therapeutic studies treatment conditions were prepared in at least triplicate within each batch including internal controls with sham manipulations. Therapeutic response was assessed via the CellTiter 96® AQueous One Solution Cell Proliferation Assay (Promega; Madison, WI) at 490nm absorbance in a BioTek® Epoch Microplate Spectrophotometer.

III Coarse-grained Model

Motivated by the preceding observations and by previous work Pisco et al. (2013); Zhou et al. (2014), we now consider a simple coarse-grained model (Fig. 2) for phenotypic heterogeneity in tumor cells. We consider that the population of cancer cells consists of two distinct subpopulations; drug-sensitive or drug-tolerant. Based on our experimental results, we denote the drug-sensitive population by $E$ (for epithelial phenotype) and the drug-tolerant population by $M$ (for mesenchymal phenotype). The processes that control the evolution of tumor heterogeneity are as follows: 1) birth: each $E$-type or $M$-type cell gives rise to birth of new cells of the same type with rates $k_{E}$ and $k_{M}$, respectively; 2) death: each $E$-type ($M$-type) cell degrades with rates $\mu_{E}$ ($\mu_{M}$); 3) phenotypic switching: an $E$ cell can switch to a $M$ cell with rate $k_{EM}$, and $M$ cell can switch back to $E$ cell with rate $k_{ME}.$

At any time $t$, the state of the system is defined by the number of $E$ and $M$ cells. The temporal evolution of the corresponding probability distribution is given by the master equation:

[TABLE]

where $P(E,M,t)$ denotes the probability that there are $E$ and $M$ numbers of epithelial and messenchymal cells present at time $t$.

Within the framework of this model (Fig. 2), we now address a key issue: How to characterize the initial heterogeneity (i.e. prior to the start of drug exposure) in tumor cells. To quantify this heterogeneity, let us consider the fraction of $M$- cells in the population. Consider the case that we isolate different sample populations from the tumor each corresponding to a fixed number ($N_{0}$) of cells, see Fig. 3. Let $p_{0}=M/N_{0}$ denote the fraction of $M$-cells in the sample population. We propose to quantify tumor heterogeneity by considering $p_{0}$ as a random variable drawn from a distribution $\rho(p_{0})$, characterized by its mean $\langle p_{0}\rangle$ and variance $\sigma^{2}_{p_{0}}=\langle p_{0}^{2}\rangle-\langle p_{0}\rangle^{2}$. Thus, initial tumor heterogeneity is characterized not just by the presence of drug-tolerant $M$ cells in the sample but also by variations in the number of $M$ cells from sample to sample, characterized by the distribution $\rho(p_{0})$.

IV Analytical results for parameter estimation

We now consider the stochastic process governing evolution of the tumor population upon treatment with drugs. Upon exposure to drugs, it is a reasonable assumption that growth is inhibited, so accordingly we set $k_{E}=k_{M}=0$. As discussed in the preceding section, we consider the evolution of different sample populations, each of which has a fixed initial size $N_{0}$ such that the fraction of $M$ cells is drawn from a distribution $\rho(p_{0})$. In this limit, the key parameters of the model are: $k_{EM},k_{ME},\mu_{E},\mu_{M},\langle p_{0}\rangle,\sigma^{2}_{p_{0}}$. In what follows, we derive analytical results that can be used to estimate model parameters by analyzing the distribution of surviving cells upon drug exposure.

We note that, within our model, the evolution of each cell in the population is independent of the state of the remaining cells. Correspondingly, we first focus on the time evolution of a single tumor cell, which is initially either $E$-type or $M$-type with corresponding probabilities as $1-p_{0}$ and $p_{0}$, Let us denote by $P_{E}(P_{M})$ the probability that the cell is $E$($M$)-type at time $t$, conditional on the initial probability $p_{0}$ for it to be $M$-type . The corresponding probability generating function for the single cell, conditional on the value of $p_{0}$ ($g(z_{1},z_{2},t|p_{0})$ $=\sum_{\eta_{E}}\sum_{\eta_{M}}z_{1}^{\eta_{E}}z_{2}^{\eta_{M}}P(\eta_{E},\eta_{M},t|p_{0})$), can be expressed as

[TABLE]

It is straightforward to derive analytic expressions for $P_{E}(P_{M})$ and to thereby obtain an expression for $g(z_{1},z_{2},t|p_{0})$ (Supplementary Material A). Now, let $G(z_{1},z_{2},t)$ denote the probability generating function corresponding to $P(E,M,t)$, the probability that we have ${E}$ and ${M}$ number of sensitive ($E$-type) and resistant ($M$-type) cells in the entire population at time $t$. Since each cell in the population (initial size $N_{0}$) evolves independently, the probability generating function for the joint distribution at time $t$ (averaging over the initial choice of $p_{0}$) is given by

[TABLE]

where $\rho(p_{0})$ is the probability distribution function for the initial fraction of $M$-type cells ($p_{0}$).

The expression derived for the generating function, Eq. (3), can be used to derive analytic expressions for all the moments of the marginal distributions corresponding to $E$-type and $M$-type cells at time $t$. For example, expressions for mean number of $E$ and $M$ cells can be obtained using $\langle E\rangle=dG/dz_{1}|_{z_{1}=1,z_{2}=1}\text{and}\langle M\rangle=dG/dz_{2}|_{z_{1}=1,z_{2}=1}$, respectively (see Supplementary Material A). This leads to the following expression for mean number of surviving cells at time $t$, $\langle N\rangle=\langle E+M\rangle$:

[TABLE]

where

[TABLE]

It is clear from the above expression that by fitting the curve corresponding to the mean number of surviving cells as a function of time, the three parameter combinations: $\alpha_{0}$, $\gamma_{0}$, and $\mu_{E}-\mu_{E}p_{0}+\mu_{M}p_{0}$ can be determined.

To extract the remaining model parameters based on time-course data, we have to turn to analytic results for the higher moments. For example, we can use expressions for the Fano factor ($F$) associated with total number of surviving cells, which is given by $F=\sigma^{2}_{N}/\langle N\rangle$ with $\sigma^{2}_{N}=\langle N^{2}\rangle-\langle N\rangle^{2}$ denoting the variance in the number of surviving cells. The expression for $\sigma^{2}_{N}$ can be obtained using

[TABLE]

where $\sigma^{2}_{E}=\langle E^{2}\rangle-\langle E\rangle^{2}$ and $\sigma^{2}_{M}=\langle M^{2}\rangle-\langle M\rangle^{2}$ are variances associated with the marginal distributions for the $E$ and $M$ cells respectively, and $C_{EM}=\langle EM\rangle-\langle E\rangle\langle M\rangle$ is the correlation between numbers of $E$ and $M$ cells. We obtain an explicit expression for the Fano factor given by (see Supplementary Material A):

[TABLE]

The Fano factor $F$ is a measure of deviations from the Poisson distribution, for which $F=1$. If $F<1$ or $F>1$, the distribution is sub-Poissonian or super-Poissonian, respectively. Before turning our attention to approaches for parameter estimation, let us first examine the expression derived for the Fano factor. We note that in the absence of initial variability in the fraction of $M$ cells (i.e. $\sigma^{2}_{p_{0}}=0$), the Fano factor of surviving population is simply given by $F=1-\langle N\rangle/N_{0}$. As the cells are not dividing due to the exposure to drugs, the mean number of surviving cells ($\langle N\rangle$) is less than the initial population of cells ($N_{0}$) for $t>0$. Thus, in this case, the Fano factor is always less than one and the distribution of cell population follows a sub-Poissonian distribution. However, given variability in the fraction of $M$-cells in the initial population, the Fano factor can potentially exceed one making the distribution super-Poissonian. This result implies that the observation of a Fano factor in excess of 1 in the distribution of surviving cells is an indicator of variance in the fraction of $M$-cells in the initial population. Thus the measurements of the moments of surviving cell populations can provide evidence for phenotypic heterogeneity in tumor populations prior to drug treatment.

To gain more quantitative insight into the initial heterogeneity, we need to estimate the parameters characterizing the mean and variance of $\rho(p_{0})$. Let us rewrite Eq. (7) in a more compact form by regrouping terms in the expression to yield the following form

[TABLE]

with

[TABLE]

Using the expressions for the mean and Fano factor of the surviving population as functions of time, Eqs. (4) and (8), we can estimate four of the parameter combinations, namely, $\alpha_{0}$, $\gamma_{0}$, $\mu_{E}-\mu_{E}p_{0}+\mu_{M}p_{0}$, and $(\sigma_{p_{0}}(\mu_{E}-\mu_{M}))/\alpha_{0}$. Correspondingly, we need additional experiments to determine the entire set of 6 model parameters. As we now show, a set of measurements that accomplish this can be obtained by starting from different initial conditions.

The proposed protocol is motivated by that fact that experimental techniques such as fluorescence-activated cell sorting (FACS) can be used to prepare the samples in specified initial states. With this in mind, we begin from an initial condition where all cells are $E$-type i.e. $p_{0}=0,\sigma^{2}_{p_{0}}=0$. Using the derived results, we can determine the parameters $\mu_{E}$, $\alpha_{0}$ and $\gamma_{0}$. We next consider the initial condition to be all $M$-type cells i.e. $p_{0}=1,\sigma^{2}_{p_{0}}=0$. Analysis of the corresponding time-course measurements of the number of surviving cells can now be used to estimate the parameter $\mu_{M}$. Having estimated values of $\alpha_{0}$, $\gamma_{0}$, $\mu_{E}$, and $\mu_{M}$, we can now find the switching rates, $k_{EM}$ and $k_{ME}$, from Eq.(5). Finally, using the expressions for the mean number of surviving cells and corresponding Fano factors for arbitrary $p_{0}$, we can get explicit expressions for the probability $p_{0}$ and variance $\sigma^{2}_{p_{0}}$ (Supplementary Material A):

[TABLE]

The above results are expressed in terms of experimentally measurable quantities, involving mean values $\langle N\rangle_{0}$ (for $p_{0}=0$) and $\langle N\rangle_{1}$ (for $p_{0}=1$) and the Fano-factor of the total surviving population, and thus can be used to estimate the population heterogeneity at the start of drug exposure based on time-course measurements of the surviving population size.

V Modeling generation of tumor heterogeneity

The analysis in the preceding section holds regardless of the source of initial heterogeneity in tumor populations. In this section, we explore how the model introduced for tumor cell dynamics can also be used to analyze a potential mechanism for generation of tumor heterogeneity. We note that the proposed model in Fig. 2 can be seen as a generalized version of the celebrated Luria-Delbrück (LD) model with the important addition that, in the present case, the transition between the two phenotypes is reversible (as opposed to the Luria-Delbrück case). However, while the LD model can be solved exactly Antal and Krapivsky (2011), the exact analytical solution of the reversible model in Fig. 2 is not known, to the best of our knowledge. Nevertheless, as we show below, exact expressions for the mean and variance of the number of $E$-type and $M$-type cells can be obtained and used to characterize heterogeneity in tumor cell populations. We can use the master equation, Eq.(1), to derive expressions for the mean number of $E$ and $M$ cells at any time $t$ (Supplementary Material B). Using these expressions, the mean number of surviving cells $\langle N\rangle=$ $\langle E\rangle+\langle M\rangle$ is given by:

[TABLE]

where

[TABLE]

with $k_{E}^{f}=k_{E}-\mu_{E}$ and $k_{M}^{f}=k_{M}-\mu_{M}$ representing the effective birth rates for $E$-type and $M$-type cells respectively, while $E_{0}$ and $M_{0}$ are the initial numbers of $E$ and $M$ cells at $t=0$.

The results show that the mean number of surviving cells at any time $t$ is characterized by six parameters: initial number of $E$ and $M$ cells ($E_{0},M_{0}$), two effective birth rates ($k_{E}^{f}$, $k_{M}^{f}$) and two switching rates ($k_{EM},k_{ME}$). Given that the initial population can be chosen in a controlled manner, we can use the results for the mean population size to determine some of the model parameters. Specifically, we can set $M_{0}=0$ as the initial condition and fitting the data to obtain the coefficient of exponential terms in Eq. (11) will yield $k_{E}^{f}$. Next, we can set $E_{0}=0$ and Eq. (11) will allow us to extract $k_{M}^{f}$. Once we estimate $k_{E}^{f}$ and $k_{M}^{f}$, we can extract the switching rates using the estimated values of $\gamma$ and $\alpha$, using Eq. (V). That is, the proposed procedure allows us to estimate the parameter combinations, $k_{E}^{f}$ and $k_{M}^{f}$ as well as the parameters $k_{EM}$ and $k_{ME}$.

In order to estimate the remaining model parameters, we need to consider the higher moments. While obtaining analytical expressions for the full probability distribution is still an open problem, higher moments can be calculated in a straightforward manner. For example, using Eq. (1) the evolution equation for $\langle E^{2}\rangle=\sum E^{2}P(E,M,t)$, $\langle M^{2}\rangle=\sum M^{2}P(E,M,t)$ and $\langle EM\rangle=\sum EMP(E,M,t)$ is given by

[TABLE]

The above set of equations can be solved to get explicit expressions for $\langle E^{2}\rangle$, $\langle M^{2}\rangle$ and $\langle EM\rangle$ at any time $t$ (Supplementary Material B), which can be used to get the variance ($\sigma^{2}_{N}$) in the total number of surviving cells by using Eq. (6). The analytic expression for the variance, in combination with the expression for mean number of surviving cells, can be used to extract all model parameters. Furthermore, the extracted parameters can then be compared with the parameters derived based on tumor cell dynamics after exposure to drugs. The comparisons can provide insight into the relative roles of adaptation and selection in driving tumor heterogeneity. The scenario wherein the switching parameters $k_{EM}$ and $k_{ME}$ are effectively unchanged upon exposure to drugs, whereas $\mu_{E}$ and $\mu_{M}$ increase favors selection as the dominant driver of tumor heterogeneity. However, significant changes in the switching rates $k_{EM}$ and $k_{ME}$ are indicative of a role for adaptation as well in the generation of tumor heterogeneity.

VI Characterizing the distribution of the fraction of resistant cells

The results derived for the moments can also be used to characterize the probability distribution $\rho(p_{0})$ for the fraction of $M$-type cells. Recall that we must have $0\leq p_{0}\leq 1$ and furthermore the first two moments of $\rho(p_{0})$ can be obtained using the procedure outlined in the previous section. Thus a natural choice to characterize the distribution $\rho(p_{0})$ is to take it to be the Beta distribution with the mean and variance as determined by measurements. Such a distribution is expressed in terms of two exponents ($\alpha$ and $\beta$) as

[TABLE]

with mean $\langle p_{0}\rangle=\frac{\alpha}{\alpha+\beta}\text{ and variance}~{}~{}\sigma^{2}_{p_{0}}=\frac{\alpha\beta}{(\alpha+\beta)^{2}(\alpha+\beta+1)}.$ The two parameters of the Beta distribution can be estimated using the experimentally determined mean and variance. Explicitly, these are given by

[TABLE]

To test this approach for characterizing the initial heterogeneity, we compare the Beta distribution with the results obtained from stochastic simulations of the model. Specifically, we carried out stochastic simulations using the Gillespie algorithm Gillespie (1977) for the model in Fig.2 starting with 200 sensitive E-type cells and no resistant M-type cells. The empirically determined distribution for the fraction of $M$-type cells ($\rho(p_{0})$) is then compared to the Beta distribution with the same mean and variance as the empirical distribution. The results obtained are shown in Fig. 4, which indicate that the Beta distribution is an excellent approximation for the range of parameters explored.

VII Discussion

To summarize, we have studied a coarse-grained stochastic model to quantify phenotypic heterogeneity in a population of cancer cells. Motivated by the experimental observation that both chemoresistance and TGF-$\beta$ induced EMT lead to similar outcomes, the model assumes that a cell has two phenotypes corresponding to whether it is drug-sensitive or drug-resistant. Importantly, the model is also consistent with epigenetic mechanisms for generating phenotypic heterogeneity in cancer, given that it allows reversible phenotypic switching between sensitive and resistant cells.

For the model considered, we have derived analytic results, both in the presence and absence of chemotherapeutic agents, which provide insights into the role of phenotypic switching in generating population heterogeneity. One of the issues that we address through these results focuses on quantifying initial heterogeneity in the fraction of resistant cells, characterized by mean fraction of resistant cells $\langle p_{0}\rangle$ in a tumor and its variance $\sigma^{2}_{p_{0}}$. For this, we propose a protocol that can be used to estimate the model parameters based on measurements of mean and variance of the surviving population of tumor cells. Furthermore, our analysis also leads to a condition, in terms of experimentally measurable quantities, whose value serves as an indicator for the presence of initial heterogeneity in the fraction of resistant cells.

While the proposed method allows us to estimate the mean and variance of the fraction of resistant cells prior to therapy, obtaining an exact analytical form for the entire distribution appears to be challenging. However, our simulation results suggest that this distribution is well approximated by the Beta distribution, which can be characterized by using the mean and variance of the surviving population. Besides characterizing initial heterogeneity in the cancer cell population, the estimated model parameters can also be useful in analyzing the complex roles of adaptation and selection in the acquisition of chemoresistance. Furthermore, the results obtained provide exact analytical expressions characterizing the distribution of of tumor cell population under treatment by drugs. Going forward we envision further model development in dialog with experiments that longitudinally monitor phenotypic changes in time lapse microscopy studies, either by quantitative analysis of morphometric parameters or implementing fluorescent reporters of EMT which have been recently developed Toneff et al. (2016). These results can serve as important inputs to future work focusing on evaluation of the hypothesis that model-informed design of treatment schedule and dose parameters may reduce the emergence of chemoresistance.

**Acknowledgments: **The authors gratefully acknowledge funding support from the NIH through grants 3U54CA156734-05S3 (as part of the UMass Boston/Dana Farber-Harvard Cancer Center U54 partnership) and from the National Cancer Institute R00CA155045 ( PI: JPC), and a Sanofi Genzyme doctoral fellowship which supported Gwendolyn Cramer. We would also like to acknowledge funding support from the UMass Boston Healey award.

Author contributions: NK, BS, JC and RVK designed research; NK, SZD and RK carried out theoretical/ computational aspects of research, GC carried out experimental aspects of research; NK, BS, JC and RVK wrote the paper and all authors reviewed the paper.

Supplementary Material A: Analytical results for surviving population upon exposure to drugs

Derivation of single cell generating function

Consider a single cell that can exist either as an epithelial($E$) or mesenchymal($M$) cell. The rates of switching between these two phenotypes ($E$ and $M$) are given by $k_{EM}$ and $k_{ME}$ and the rates of cell death are given by $\mu_{E}$ and $\mu_{M}$, respectively. Here we assume no cell divisions (i.e no new production of cells) for either $E$ or $M$ due to high levels of external drugs. For such a system, we consider first the temporal evolution of a single cell, given the initial probability $p_{0}$ of the cell being $M$-type. The corresponding probability generating function is $g(z_{1},z_{2},t|p_{0})$ $=\sum_{\eta_{E}}\sum_{\eta_{M}}z_{1}^{\eta_{E}}z_{2}^{\eta_{M}}P(\eta_{E},\eta_{M},t|p_{0})$, where $\eta_{E}$ and $\eta_{M}$ can have values 0 or 1, and $P(\eta_{E},\eta_{M},t|p_{0})$ is the probability of having $\eta_{E}$ and $\eta_{M}$ number of cells at time $t$, given the initial proability $p_{0}$ of the cell being $M$-type. Clearly $P(1,1,t|p_{0})=0$ since we are starting with a single cell and no new cells are created. Correspondingly at any time $t$, we have only three possibilities, either $\eta_{E}=1$ and $\eta_{M}=0$ or $\eta_{E}=0$ and $\eta_{M}=1$ or $\eta_{E}=\eta_{M}=0$. Thus we obtain

[TABLE]

Denote $P(1,0,t|p_{0})=P_{E}(t)$, $P(0,1,t|p_{0})=P_{M}(t)$ and $P(0,0,t|p_{0})=P_{0}(t)$, and then using the normalization condition, $P_{E}(t)+P_{M}(t)+P_{0}(t)=1$, we obtain the single particle generating function as derived in the main text, Eq. (2).

To find an explicit expression for the single cell generating function $g(z_{1},z_{2},t|p_{0})$, we need to find expressions for the probabilities $P_{E}(t)$ and $P_{M}(t)$. For this, we begin from their evolution equations:

[TABLE]

These equations can be solved to give the following expressions for the temporal evolution of the probabilities:

[TABLE]

where $\alpha_{0}$ and $\gamma_{0}$ denote combinations of the model parameters, as given in Eq.(5) in the main text and $p_{0}=P_{M}(t=0)$ as discussed above. Using the expressions for $P_{E}$ and $P_{M}$ from Eq. (Derivation of single cell generating function) in Eq. (2), we obtain the single particle probability generating function.

Derivation of mean and Fano factor for surviving population

Here we provide the details for the derivation of expressions for temporal evolution of moments associated with total surviving population. For this, we start from the expression for the generating function,

[TABLE]

and using Eqs.(A4), (2), (Derivation of single cell generating function) and denoting $\langle p_{0}\rangle=\int p_{0}\rho(p_{0})dp_{0}$ as the mean value of $p_{0}$, arrive at the following expressions for the mean numbers of epithelial and mesenchymal cells,

[TABLE]

where

[TABLE]

and $\alpha_{0}$ and $\gamma_{0}$ are given by Eq. (5) in the main text. Using Eq.(Derivation of mean and Fano factor for surviving population), we see that the mean number of total surviving cells ($E+M$) is given by

[TABLE]

We turn next to the second moment and, using the same generating function, the variances in $E$ and $M$ cells are given by

[TABLE]

leading to the following expression for the Fano factor:

[TABLE]

Substituting the expressions for $S_{1}$ and $Q_{1}$ and simplifying these Fano factors can be reexpressed as

[TABLE]

Note that both the Fano factors, $F_{E}$ and $F_{M}$, are always less than one if $\sigma^{2}_{p_{0}}=0$. That is, a Fano factor greater than 1 is an indication of the presence of initial variability in the fraction of $M$ cells in the population.

To derive an expression for the variance of the total population $N=E+M$, we use the relation $\sigma^{2}_{N}=\sigma^{2}_{E}+\sigma^{2}_{M}+2C_{EM}$, where $C_{EM}=\langle EM\rangle-\langle E\rangle\langle M\rangle$ is the correlation between $E$ and $M$, and $\langle EM\rangle$, in terms of the generating function, is

[TABLE]

The expression for the Fano factor for the total population $N$ can be written as

[TABLE]

which, on using Eq. (Derivation of mean and Fano factor for surviving population), can be rewritten as Eq. (7) in the main text.

Expressions for $p_{0}$ and $\sigma^{2}_{p_{0}}$

As discussed in the main text, we propose a three-step procedure. First, we set $p_{0}=0$, which, using Eqs. (A7) and (A12), gives

[TABLE]

and next, setting $p_{0}=1$,

[TABLE]

Finally, using the expressions for mean number of surviving cells, Eq.(A7), and corresponding expression for the Fano factor, Eq. (A12), for arbitrary $p_{0}$, we get explicit expressions for the probability $p_{0}$ and variance in initial M-type cells $\sigma^{2}_{p_{0}}$. The resulting expressions are in terms of mean values $\langle N\rangle_{0}$ (for $p_{0}=0$) and $\langle N\rangle_{1}$ (for $p_{0}=1$) and the mean and Fano factor for a given arbitrary $p_{0}$, as shown in the main text.

Supplementary Material B: Analytical results in the growth phase

In this section, we provide details of the derivation of the moments for the surviving cell populations in the growth phase.

First moments

Multiplying Eq. (1) by E or M and summing over all possible values of E and M, we obtain the evolution equations for $\langle E\rangle=\sum EP(E,M,t)$ and $\langle M\rangle=\sum MP(E,M,t)$ :

[TABLE]

These equations can be solved to get:

[TABLE]

where $E_{0}$ and $M_{0}$ are initial values for the number of $E$-type and $M$-type cells respectively, and

[TABLE]

with

[TABLE]

representing the effective birth rates for $E$-type and $M$-type cells, respectively. For the total population, $N=E+M$, using Eq. (First moments), the mean $\langle N\rangle=$ $\langle E\rangle+\langle M\rangle$ can then be written as:

[TABLE]

Second moments

Using Eq. (1), we can write the evolution equation for

[TABLE]

as:

[TABLE]

Eliminating variables in Eq. (Second moments), we obtain a single ODE for $\langle E^{2}\rangle$,

[TABLE]

where $c_{0}$, $c_{1}$ and $c_{2}$ are functions of the model parameters

[TABLE]

and $\mathcal{K}(t)$ is a known function in $t$:

[TABLE]

with

[TABLE]

The solution for $\langle E^{2}\rangle$ in Eq. (B7) is:

[TABLE]

where $\lambda_{i}$, $i=1,2,3$ are the roots of

[TABLE]

and $s_{1},s_{2},s_{3}$ are constants to be determined from the initial condition. Once we have obtained an expression for $\langle E^{2}\rangle$, we can readily derive expressions for $\langle M^{2}\rangle$ and $\langle ME\rangle$:

[TABLE]

Taking suitable initial conditions, we can find the constants $s_{1},s_{2},s_{3}$ that will help specify the temporal evolution of the second moments in terms of the model parameters.