I. INTRODUCTION

Advances in single-cell technologies have exposed remarkable differences in phenotype and expression patterns between individual cells within the same isogenic cell population [1]–[9]. While some of this variation can be linked to extrinsic factors (i.e., cell-cycle stage, cell size, local extracellular environment), increasing evidence points to the role of stochastic processes inherent to gene expression in driving random fluctuations (noise) in gene product levels [10]– [24]. This intercellular phenotypic heterogeneity can play important functional roles in diverse biological processes, from driving genetically-identical cells to different cell fates [25]–[35] to allowing microbes and cancer cells to hedge their bets against uncertain environmental changes [36]–[47].

While single-cell sequencing tools can probe phenotypic heterogeneity within a given cell population, they only provide a static picture of different cell states. Characterizing the dynamics of individual cells transitioning between of cellG-roswctolaonietse

from single cells Grow colonies from single cells different states with multi-generational time scales remains a fundamental challenge in advancing the field of singlecell biology. The Luria-Delbrück experiment, also called the “Fluctuation Test", demonstrated that genetic mutations arise randomly in the absence of selection – rather than in response to selection – and led to a Nobel Prize [48]. We leverage a similar fluctuation-style assay to infer switching dynamics between cellular states. The proposed methodology relies on first growing single cells into colonies and then measuring each colony’s cell-state composition (i.e, the number of cells in different cell states). As a simple example, drug treatment of a single-cell derived colony of cancer cells can classify individual cells into two phenotypic states: drugsensitive (drug-tolerant) cells that die (survive) in response to treatment. Repeating this process for several colonies yields a distribution for the number of surviving cells. How can we exploit this measured distribution to understand the interchange between cell states?

Starting with a system of reversible switching between two cell states, we derive three different analytical formulas employing different approximations to quantify fluctuations in cell-state composition across colonies and benchmark them with exact stochastic simulations of the cell proliferation/switching process. From a practical standpoint, these formulas are a valuable tool to back-calculate the switching rates from experimentally measured fluctuations at a single time point of colony expansion. In the context of the cancer example, inferring the timescale of switching between drugsensitive and drug-tolerant states is critical for the design of optimal drug treatment schedules [49]. Later on, we further generalize these results to multiple cell states with arbitrary switching topologies. We start by reviewing the original Luria-Delbrück experiment done 80 years ago.

Cells switch between 2 states

By the early 20th century it was known that bacteria can acquire resistance to infection by phages (bacterial viruses). However, it was debated whether mutations leading to resistance were directly induced by the virus (Lamarckian theory), or if they arose randomly in the population before viral infection (Darwinian theory). To discriminate between these alternative hypotheses, Luria & Delbrück designed an elegant experiment where single cells were isolated and grown into colonies. After allowing the colonies to grow for some duration, they were infected by the T1 phage, and the number of resistant bacteria were counted across colonies. If mutations are virus-induced (i.e., no genetic heritable component to resistance), then each bacterium has a small and independent probability of acquiring resistance, and the colony-to-colony fluctuations in the number of resistant cells should follow Poisson statistics. In contrast, if mutations occur randomly before viral exposure, then the number of surviving bacteria will vary considerably across colonies depending on when the mutation arose in the colony expansion (Fig. 1). The data clearly showed a non-Poissonian skewed distribution for the number of resistant bacteria, validating the Darwinian theory of evolution [48].

The Luria-Delbrück experiment that came to be known as the “Fluctuation Test", not only addressed a fundamental evolutionary question leading to a Nobel Prize, but also laid the foundations for the field of bacterial genetics. Apart from its biological significance, the fluctuation test exemplifies the usage of stochastic analysis for uncovering hidden processes even though the underlying cell states may not be directly observable. Publication of the LuriaDelbrück article in 1943 catalyzed rich theoretical work deriving probability distributions for the number of resistant cells based on different biological assumptions [50]–[52], and led to statistical methods for estimating mutation rates from fluctuations in the data [53]–[55]. We refer the reader to [56] for an excellent review of mathematical developments related to the Luria-Delbrück experiment.

The fluctuation test was recently used to study cancer drug resistance [57]. Individual melanoma cells were isolated from a clonal cell population by single-cell FACS sorting, and then grown into colonies. After allowing single cells to expand for a few weeks, the colonies were treated with a targeted cancer drug, vemurafenib. Intriguingly, the colonyto-colony fluctuations in the number of surviving cells were significantly larger than a Poisson distribution, but an order of magnitude smaller than what is predicted by the mutation State 1

State 2 Fast-switching

Medium-switching model. Subsequent analysis showed that stochastic expression of several resistant markers drives individual cells to rNevoenr-sreibsilsytanswtitch between drug-sensitive and drug-tolerant sRtaetseisstapnrtior to drug exposure [58], and the latter stateGrcoawncolonies transform into a stably-resistant state in the presence ofrfotmhesingle cells drug [57]. We next present the mathematical framework for modeling cell-state transitions in an expanding cell colony, followed by the derivation of formulas quantifying the exteAndt d T1 Phage of fluctuations as a function of the switching rates.

III. FLUCTUATION-APPROACH TO INFER CELL-STATE

SWITCHING

Consider a scenario as in Fig. 2 where cells within a population can reside in two states (State 1 & 2). Cells proliferate and reversibly switch between states, and the rates of switching determine the transient heritability of a state, i.e., IhnodwucemdManuytatgioennerations it takes tSoposnwtaintcehoubsaMcuktattoionthe other state. Let fˆ2 denote the average fraction of cells in State 2 in the original population. Single cells are randomly drawn from the population (through serial dilutions or FACS sorting), and expanded into colonies. Note that the state of the starting single cell is unobservable, as we only consider a single endpoint measurement. After growing the colonies for a certain duration of time, each colony is assayed for the fraction of cells in State 2 (or State 1). The basic idea is that if switching between states is relatively fast (several switches happen in the growth duration), then the fraction of State 2 cells will rapidly equilibrate to fˆ2 in each colony, and colony-to-colony fluctuations will be minimal (Fig. 2). In contrast, if switching is slow, then based on the memory of the initial cell, colonies will primarily be composed of cells in either State 1 or 2 generating large colony-to-colony fluctuations (Fig. 2). In essence, fluctuations in colony cellstate composition reveals the timescale of switching, with slower relaxation kinetics driving higher fluctuations.

To directly relate these fluctuations to the switching kinetics, we model the cell-colony expansion process by considering that cells proliferate with a constant rate kx implying a mean cell-cycle time of 1/kx (i.e., time taken by each cell to finish cell-cycle and make two daughters). Cells in State 1 transition to State 2 with a constant rate k1, and switch back to State 1 with a constant rate k2 resulting in the average fraction fˆ2 := k1 . ( 1 )

k1 + k2 We make several simplifying assumptions in our model formulation: • The proliferation rate of a cell is the same irrespective of the cellular state. • There is no cell death in the expanding colony. • The population remains in the exponential growth phase during the time span of the experiment. • The switching rates are independent of the number of cells in the colony.

Assuming that the starting cell at time t = 0 is chosen randomly from the original population, its state follows a Bernoulli random variable: the cell is in State 2 with probability fˆ2 and in State 1 with probability 1 − fˆ2. Let stochastic processes x1(t) and x2(t) denote the number of cells in State 1 and 2, respectively, with total cell number x(t) = x1(t) + x2(t). Our goal is to quantify colony-to-colony fluctuations in the fraction of cells in State 2 f 2(t) = x2(t) , ( 2 )

x(t) at time t of colony expansion. It is important to point out that we work with the fraction of cells in State 2, and not the number of cells in State 2. The rationale for this is that fluctuations in the former appear more robust to the cellcycle time distribution. For example, stochastic simulations of an expanding colony with an exponentially-distributed or a lognormally-distributed cell-cycle time (with the same mean) result in different levels of fluctuations in x2(t), but similar fluctuations in f 2(t) (Fig. 3). Employing different assumptions, we next derive closed-form formulas for the coefficient of variation (standard deviation over mean) of f 2(t) across colonies.

Ignoring any stochasticity in the cell proliferation and switching processes, colony growth from a single cell can be modeled as an ordinary differential equation dx dt = kxx, dx2 = kxx2 + k1x − (k2 + k1)x2 dt =⇒ d f2 = k1(1 − f2) − k2 f2. ( 3 ) dt Solving ( 3 ) for a Bernoulli initial condition f2(0) = 1 with probability fˆ2 (starting cell is in State 2) and f2(0) = 0 with probability 1 − fˆ2 (starting cell is in State 1) yields f2(t) = fˆ2 + (1 − fˆ2)e−(k1+k2)t with probability fˆ2 f2(t) = fˆ2 − fˆ2e−(k1+k2)t

with probability 1 − fˆ2.

This results in the following mean fraction h f2(t)i = fˆ2 fˆ2 + (1 − fˆ2)e−(k1+k2)t + (1 − fˆ2) fˆ2 − fˆ2e−(k1+k2)t = fˆ2 where h i denotes the expected value of a random process. As intuitively expected, at any time in the lineage expansion the mean fraction of cells in State 2 across colonies is the same as that in the original population. Using a similar approach for the second-order moment h f22(t)i = fˆ2 fˆ2 + (1 − fˆ2)e−(k1+k2)t 2 + (1 − fˆ2) fˆ2 − fˆ2e−(k1+k2)t 2 =⇒ h f22(t)i = fˆ22 + fˆ2(1 − fˆ2)e−2(k1+k2)t .

This leads to the following colony-to-colony fluctuations

CV22(t) := h f22(t )fˆ2i2 − fˆ22 = 1 −fˆ2 fˆ2 e−2(k1+k2)t ( 8 ) with CV2 being the coefficient of variation (CV ) of the State 2 fraction. The formula shows that for a fixed fˆ2 and time t, faster switching attenuates fluctuations with CV2 → 0 as k1, k2 → ∞. By defining a dimensionless quantity we can rewrite ( 8 ) as

Z := 1 −fˆ2 fˆ2 e− 2kZxt .

( 4 ) ( 5 ) ( 6 ) ( 7 ) ( 9 ) Given a time t of CV2 measurement, switching rates can be inferred by simultaneously solving ( 1 ) and ( 10 ) based on apriori knowledge of the cell proliferation rate, and fˆ2 is the average fraction of State 2 cells across colonies.

To capture low-copy number effects at early stages of colony expansion, we consider a stochastic model formulation where now x(t), x1(t), x2(t) ∈ {0, 1, 2, . . .} are integervalued random processes. For defining the model and its subsequent analysis we work with x(t) and x2(t). The time evolution of these populations is governed by the events in Table 1 that occur probabilistically with rates given in the third column. Whenever the event occurs, the corresponding reset map is activated and cell numbers increase/decrease by one. This formulation corresponds to the cell-cycle time being drawn from an exponential distribution with mean 1/kx. Moreover, the time spent in States 1 and 2 is exponentiallydistributed with means 1/k1 and 1/k2, respectively.

Given that the event rates are linear functions of the statespace, the statistical moments of x(t), x2(t) can be obtained exactly using standard tools from moment dynamics. In particular, the time derivative of the expected value of any continuously differentiable function ϕ(x, x2) is given by dhϕ(x, x2)i

= hg(x, x2)i, dt g(x, x2) = kx(x − x2) (ϕ(x + 1, x2) − ϕ(x, x2)) + kxx2 (ϕ(x + 1, x2 + 1) − ϕ(x, x2)) + k1(x − x2) (ϕ(x, x2 + 1) − ϕ(x, x2)) + k2x2 (ϕ(x, x2 − 1) − ϕ(x, x2)) Substituting the solution of ( 12 )-( 13 ) in ( 16 ) results in h f2(t)i = fˆ2 and the following formula for the extent of fluctuations (16a) (16b) ( 18 ) ( 19 ) Based on the Bernoulli cell-state assignment of the starting cell, solving the linear dynamics system ( 12 )-( 13 ) with initial conditions hx(0)i = 1, hx2(0)i = fˆ2, determines the statistical moments of x(t) and x2(t) at any time in the colony expansion. How do we use the moments of x(t) and x2(t) to obtain moments of their ratio f 2(t) = x2(t)/x(t)? We next discuss two alternative approaches for doing this.

Analyzing the moment dynamic equations reveals hxx2i − hxihx2i = 1 − e−kxt > 0 ( 15 )

hxihx2i a positive covariance between the number of cells in State 2 and the colony size. Equation ( 15 ) points to weakly correlated x and x2 = x f2 for short times, implying an even weaker correlation between x and f2. This motivates an approximation where the fraction of cells in State 2 is independent of the colony size. Working along these lines, assuming f2(t) and x(t) to be independent allows an analytical derivation for the extent of fluctuations in f 2(t) in spite of the nonlinear dependence. To see this, recall that by definition x2 = f 2x, assuming h f2xi ≈ h f2ihxi and h f22x2i ≈ h f22ihx2i, the moment of f2 can be obtained as hx2i = h f2xi ≈ h f2ihxi =⇒ h f2i ≈ hx22i = h f22x2i ≈ h f22ihx2i =⇒ h f22i ≈ hhxx222ii . [59]–[62]. By setting ϕ(x, x2) = x and x2 in ( 11 ), one obtains the average population dynamics dhxi dhx2i

= kxhxi, = kxhx2i + k1hxi − (k1 + k2)hx2i ( 12 ) dt dt that is identical to the deterministically-formulated dynamics ( 3 ). Following the same steps for ϕ(x, x2) = x2, xx2 and x2 2 dhx2i

dt dhx22i

dt dhxx2i dt = kxhxi + 2kxhx2i = k1hxi − 2k2hx22i + k2hx2i + 2k1hxx2i − 2k1hx22i − k1hx2i + 2kxhx22i + kxhx2i = −k2hxx2i + k1hx2i − k1hxx2i + 2kxhxx2i + kxhx2i. ( 11 ) (13a) (13b) (13c) and if one grows the colony long enough, then limt→∞ CV22 = 0 as f 2(t) converges to fˆ2 in each colony. For a fixed fˆ2 and t, CV2 monotonically increases with increasing time spent in State 2, i.e., slower switching results in higher fluctuations in the fraction of State 2 cells (Fig. 4).

Here the dimensionless constants tkx and kx/k2 have important biological interpretations: tkx represents the average number of generations of colony expansion, and kx/k2 is the average number generations a cell remains in State 2 before switching to State 1. As expected from the Bernoulli initial condition, limCV22 = t→0 1 − fˆ2 fˆ2 , Moreover, from ( 20 ) h f2i = D x2 E ≈ hx2i x hxi

= fˆ2. xx2 − hhxx2ii ( 23 ) which after substituting the population number moments obtained from solving ( 12 )-( 13 ) yields the following formula CV22(t) = e−kxt (1 + 2fkˆ2xt)(1 − fˆ2) , Z = 2.

Z − 2 2Ze(− 2kZxt ) − (2 + Z)e−kxt ! 1 − fˆ2 fˆ2 In summary, we have now developed three different formulas given by ( 10 ), ( 17 ) and ( 24 ), quantifying fluctuations in fraction State 2 cells across colonies. It is interesting to note that while all formulas reduce to ( 19 ) at t = 0, they have qualitatively different initial slopes lim dCV22 = − 2kx(1 − fˆ2) < 0 (Deterministic Proliferation) t→0 dt fˆ2Z (25a) (25b) (25c) lim dCV22 = 0 (Independent Variable) t→0 dt tl→im0 dCdVt22 = kx(1fˆ−2 fˆ2)

> 0 (Small Noise) with fluctuations increasing with time in the early phase of colony expansion in (25c). All formulas share the common feature of limt→∞ CV2 = 0.

Fig. 5 shows the accuracy of these formulas by comparison with exact CV2 values as obtained by performing stochastic simulations of the system shown in Table I. While the formula based on the deterministic formulation greatly underestimates fluctuations in f2(t), the independent variable assumption that takes into account stochastic proliferation provides a better approximation, especially at initial time points. The formula based on the small noise approach goes in the wrong direction at the start and overestimates the noise, but converges to the exact value at later time points when CV2 is small.

IV. STATE TRANSITIONS BETWEEN MULTIPLE CELL

STATES

Having discussed an innovative approach to identify rates of switching between two states from measured inter-colony fluctuations, we now generalize this approach to multiple cell states. With more than two cell states there will be many different switching topologies, and one would like to identify plausible topologies from measured fluctuations in cell-state compositions. Here the fluctuation data will also be richer — one will measure both the variances and covariances in the fraction of different states across colonies. In this section, we show how to analytically predict these fluctuations for switching between an arbitrary number of cell states.

Notation: Consider a system with n ≥ 2 cell states - States 1, 2, 3 . . . , n. In each state the cell proliferates with rate kx and ki j, i 6= j is the switching rate from state i to j. We denote by fˆi as the average fraction of cells in State i ∈ {1, . . . , n}. Random processes xi(t) and fi(t) represent the number and fraction of State i cells in the expanding colony, respectively, with total cell numbers x(t) = ∑ xi(t). Our goal is to quantify the colony-to-colony fluctuations in fraction State i cells CVi2 = h fi2i − h fii2 h fii2 and the normalized covariances

Ci j = h fi f ji − h fiih f ji h fiih f ji

, i 6= j, i, j ∈ {1, 2, 3} over time as a function of the proliferation/switching rates. Having described the notation used in this section, we briefly discuss the methodology for obtaining these fluctuations using the assumptions discussed in the previous section.

Let f (t) = [ f1(t), . . . , fn(t)]T denote a column vector of the fraction of cells in each state. Starting from a single cell that is in state i with probability fˆi, modeling colony expansion as an ordinary differential equation, f (t) evolves as per the linear dynamical system d f

= A f dt where elements of the A matrix are given by

n aii = − ∑ ki j, j=1

ai j = k ji for i 6= j.

To account for stochastic proliferation one can build a model similar to Table I that enumerates all probabilistic events corresponding to cell proliferation in each state and switching between states. Let vector μ consists of all first- and second-order moments of integer-valued random processes x1(t), . . . , xn(t). Then, its time evolution can be derived as a linear dynamical system dμ

= Aμ μ. ( 31 ) dt using moment dynamic tools [59], [61]. Solving this system with initial condition hxi(0)i = fˆi, hxi2(0)i = fˆi, hxi(0)x j(0)i = 0 for i 6= j ( 32 ) ( 26 ) ( 27 ) ( 28 ) ( 29 ) ( 30 ) determines the time evolution of all the first- and secondorder moments of cell numbers in different states. Having obtained these number moments, the moments of state fractions are obtained as h fi2i ≈ hhxxi22ii , h fi f ji ≈ hxix ji for i 6= j.

hx2i ( 33 ) assuming that fi(t) is independent of x(t).

As derived in ( 23 ), for small deviations in cell numbers from their means the fluctuations in the fraction of State i is given by CVi2 ≈ hxi2i − hxii2 hxii2 + hx2i − hxi2 hxi2 − 2 (hxixi − hxiihxi) hxiihxi . ( 34 ) Using a similar Taylor series expansion approach it can also be shown that

Ci j := ≈ − h fi f ji − h fiih f ji

h fiih f ji hxix ji − hxiihx ji

hxiihx ji hxixi − hxiihxi hxiihxi + hx2i − hxi2

hxi2 − hx jxi − hx jihxi hx jihxi ( 35 ) in the limit of small fluctuations in cell numbers. Substituting moments computed from ( 31 ) in ( 34 ) and ( 35 ) provides the extent of fluctuations and covariances between state fractions over time.

We illustrate these approximations on a ring topology of state switching where k21 = k32 = k13 = 0 (Fig. 6). Comparison with exactly-computed fluctuations from stochastic simulations shows that the independent variable approximation generally provides a good approximation at earlier time points when noise levels are high, while the small noise approximation works best at later time points.