Stochastic modeling of phenotypic switching and chemoresistance in cancer cell populations

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.

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 M E and the rates of cell death are given by µ E and µ 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 ) = η E η M z η E 1 z η M 2 P (η E , η M , t|p 0 ), where η E and η M can have values 0 or 1, and P (η E , η M , t|p 0 ) is the probability of having η E and η 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 η E = 1 and η M = 0 or η E = 0 and η M = 1 or η E = η M = 0. Thus we obtain 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: These equations can be solved to give the following expressions for the temporal evolution of the probabilities: where α 0 and γ 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. (A3) 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, and using Eqs.(A4), (2), (A3) and denoting p 0 = p 0 ρ(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, where and α 0 and γ 0 are given by Eq. (5) in the main text. Using Eq.(A5), we see that the mean number of total surviving cells (E + M ) is given by We turn next to the second moment and, using the same generating function, the variances in E and M cells are given by leading to the following expression for the Fano factor: Substituting the expressions for S 1 and Q 1 and simplifying these Fano factors can be reexpressed as Note that both the Fano factors, F E and F M , are always less than one if σ 2 p0 = 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. The analytical predictions for the mean and Fano factor of the surviving E and M populations are shown in Fig. S1 for a range of parameters, along with results from stochastic simulations using the Gillespie algorithm.
To derive an expression for the variance of the total population N = E + M , we use the relation σ 2 N = σ 2 E + σ 2 M + 2C EM , where C EM = EM − E M is the correlation between E and M , and EM , in terms of the generating function, is The expression for the Fano factor for the total population N can be written as which, on using Eq. (A6), can be rewritten as Eq. (7) in the main text.
Expressions for p0 and σ 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 and next, setting p 0 = 1, 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 σ 2 p0 . The resulting expressions are in terms of mean values N 0 (for p 0 = 0) and N 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.
These equations can be solved to get: where E 0 and M 0 are initial values for the number of E-type and M -type cells respectively, and representing the effective birth rates for E-type and M -type cells, respectively. For the total population, N = E + M , using Eq. (B2), the mean N = E + M can then be written as: exp − t 2 (γ + α)