Inferring Cell-State Transition Dynamics from Lineage Trees and Endpoint Single-Cell Measurements.

Kin Correlation Analysis

(KCA) validation by comparing inferred two-state switching dynamics with direct time-lapse analysis To experimentally apply KCA, we proceed in two stages. First, we validate the method by analyzing switching between two distinct states of Esrrb expression in mouse ES cells. Second, we broaden the analysis to determine the transition dynamics of a larger set of ES cell states. Esrrb is a transcription factor central to maintaining the naïve pluripotent state, and it plays a critical role in the core pluripotency network ( Festuccia et al., 2012 ; Martello et al., 2012 ; Singer et al., 2014 ; van den Berg et al., 2008 ). Esrrb up-regulation has also been shown to facilitate fibroblast reprogramming to the induced pluripotent state ( Feng et al., 2009 ). Most importantly here, Esrrb expression in LIF+Serum culture conditions is bimodal, with cells switching between high and low expression states ( Singer et al., 2014 ). We constructed a knock-in fluorescent reporter for Esrrb expression ( Figure 2A ), and validated the reporter using smFISH (S2A–B). We acquired time-lapse movies ( Fig 2B,i ), using custom software to track individual cells over time and establish the pedigrees (lineage trees) of individual colonies (see STAR Methods). At the end of movie (~48 hours), we fixed the cells and acquired smFISH measurements of Esrrb expression ( Figure 2B,ii–iiii ). Finally, we combined the measurements, assigning smFISH Esrrb expression levels at the final time-point to the corresponding leaves of the tree ( Fig 2D,i ). Altogether, we analyzed 14 trees (299 cells) for this analysis. Consistent with previous results, Esrrb exhibited a bimodal distribution of mRNA copy number by smFISH ( Fig. 2C ) ( Kumar et al., 2014 ; Singer et al., 2014 ). To understand this distribution, we first note that a single state is expected to generate a distribution of mRNA copy numbers in individual cells, due to the stochastic, “bursty” nature of transcription and mRNA degradation, as shown previously ( Elowitz et al., 2002 ; Friedman et al., 2006 ; Ozbudak et al., 2002 ; Peccoud and Ycart, 1995 ; Suter et al., 2011 ). For many genes and cell types, the distribution of mRNA copy number is well-fit by a negative binomial distribution ( Friedman et al., 2006 ; Raj et al., 2006 ). This distribution is generated when there is a constant probability per unit time of initiating a transcriptional burst, and the number of mRNAs produced per burst follows an exponential distribution. For a gene with two expression states, we expect each state to generate a negative binomial distribution of mRNA with different burst rate and burst size parameters. These two distributions will, in general, overlap. Thus, the bimodal distribution can be explained as a linear combination of two negative binomial distributions, one for each expression state. We fit the observed Esrrb distribution to a linear combination of two negative binomial distributions. Using this fit, we assigned each cell a probability of being in either the high or low Esrrb expression state given its observed transcript count (see STAR Methods). Thus we obtained a probabilistic endpoint state assignment for each cell on each pedigree ( Figure 2D,i ). Because of the overlap between the transcript count distributions of the two Esrrb states, many cells have approximately equal probability of being in either state. The KCA framework is compatible with these probabilistic state assignments. When computing the correlation matrix ( Box 1 ), we account for probabilistic state assignments by summing over all possible pairs of states, for each pair of cells, weighting each state pair by its relative probability. When the state assignments are more ambiguous, a larger number of observations (pedigrees) is required to ensure accurate inference of transition rates (see STAR Methods for details). To infer the rates at which cells switch between Esrrb states, we analyzed these trees with KCA. We first computed pair correlation matrices for lineage distances u , ranging from 1 to 4 ( Fig. 2D,ii ). As expected, the frequency of observing two cells in the same Esrrb state decreased with increasing lineage distance. Next, using KCA we computed the switching rates that would give rise to the observed correlation matrices. For stationary Markovian dynamics, these rates should not depend on the lineage distance from which the correlation matrix is computed ( Box 1 ). The inferred rate of switching from the Esrrb low state to the Esrrb high state was 0.09±0.03 per cell cycle (errors are the standard deviation as estimated by bootstrap, see STAR Methods). The reverse rate was 0.08±0.02, per cell cycle. These rates remained constant across lineage distances of u = 1 to 4, consistent with stationary Markovian dynamics ( Fig 2D,iii–iv ). We note that the constant inferred rates imply an exponential waiting time between state transitions, a property of Markovian dynamics. To independently validate these inferred Esrrb switching dynamics, we next analyzed data from the Esrrb knock-in fluorescent reporter. We extracted the total fluorescence from each cell over the duration of the movies. Because the H2B-mCitrine is stable, its abundance diminishes only through dilution during cellular division events. These dilution events correspond to approximately halving of the fluorescent readout from each cell across cell divisions, as evident in the saw-tooth pattern of the traces shown in Fig. 2Eii–iii . We therefore focused on the “promoter activity”, or the rate of accumulation of total fluorescence (slope of the fluorescence traces shown in Fig. 2Eii–iii ), which should be proportional to the abundance of mRNA in the cell at any given time ( Singer et al., 2014 ). Indeed, the Esrrb production rate in the final cell cycle of the movie was strongly correlated with Esrrb transcript counts measured at the end of the movie, but not with that of β -actin , a homogenously expressed housekeeping gene ( Fig S2B ), providing an internal validation of both readouts. To classify Esrrb promoter activity as either high or low, we implemented a threshold on production rate at each time-point throughout the cells’ lineage history ( Fig. 2E,i ). Cell state transitions were defined as a change in the promoter activity across the threshold that persisted for at least one cell cycle (see STAR Methods). Examples of transitions can be observed in plots of total fluorescence trajectories, as shown in Fig 2E,ii–iii . Transitions from Esrrb low to high states, or high to low states, occurred with rates of 0.10±0.01 and 0.08±0.01 per cell cycle, respectively ( Fig 2E,iv ), consistent with the values inferred by KCA above (errors are the uncertainty in the observed frequencies due to finite number of observations). Finally, although this has no bearing on using KCA for inferring the transition rates, we also checked whether state transitions were more likely to occur at one particular point of the cell cycle. However, analysis revealed no strong cell cycle dependence in these data ( Fig 2F ). Together, these results suggest that the KCA method can correctly infer the reversible state switching dynamics of Esrrb , which appear to be consistent with a constant switching rate per unit time.

STAR Methods CONTACT FOR REAGENT AND RESOURCE SHARING The Lead Contact MBE is willing to distribute all materials (including constructs and engineered cell lines), datasets, software and analysis tools, and protocols used in the manuscript. Requests should be made directly to Michael B. Elowitz at melowitz@caltech.edu or by mail at California Institute of Technology. 1200 E. California Blvd., MC 114-96. Pasadena, CA 91125.

EXPERIMENTAL MODEL AND SUBJECT DETAILS Cell Line

Construction and Tissue Culture E14 cells

(E14Tg2a.4) obtained from Mutant Mouse Regional Resource Centers were used as the base line for all cell line construction. Knock-In reporters were generated using CRISPR/Cas9 with guides targeting the C-terminus of the genes of interest ( Supplementary Table 2 ), using donor vectors harboring ±300bp homology to the target locus flanking a T2A-H2B-XFP-P2A-PuroR. Single clones were first grown in 2i and isolated based on puromycin resistance and characterized for correct targeting using qPCR for genomic copy number, and then by a co-localization test of the endogenous targeted gene and XFP by smFISH ( Figs. S2A and S2C ). For the Zscan4 reporter, the 2570 base pairs upstream of the Zscan4c start codon were used as a promoter fragment reporter (as described in Zalzman, et al, 2010 ), to drive expression of H2B-mTurquoise2 on a PiggyBac integrated vector, which also contained a separate Blasticidin resistance cassette under an SV40 promtoer. Cells were maintained at 37°C and 5% CO 2 in GMEM, 10% FBS, 2 mM L-glutamine, 100 units/ml penicillin, 100 ug/ml streptomycin, 1 mM sodium pyruvate, 1000 units/ml Leukemia Inhibitory Factor (LIF, Millipore), 1X Minimum Essential Medium Non- Essential Amino Acids (MEM NEAA, Invitrogen) and 50 uM β-Mercaptoethanol. Cell lines were also stably integrated with a PiggyBac-pGK-palmitoylated-mTurquoise2/HygroR to enable 3D segmentation of cell membranes. METHOD DETAILS Time Lapse Microscopy and single-molecule Fluorescence in situ Hybridization (smFISH) Imaging For movies, cells were plated on Laminin-511 (BioLamina) in 24-well glass bottom plates (MatTek) six hours prior to the start of the movie at a density of 2000/well. Snapshots were taken at 12 minute intervals for ~48 hours, and tracked and segmented using home-grown Matlab scripts. Immediately following the end of the movie, cells were fixed in 4% Formaldehyde for five minutes at room temperature, and permeabilized in RNAse-free 70% ethanol and stored at −20°C overnight. The following day, cells were hybridized for smFISH overnight at 30°C, where genes of interest were simultaneously targeted with up to 48 20mer DNA oligos, with each gene’s probeset coupled to Alexa 555, 594, or 647 (Lifetech). Each 20mer oligo was used at ~3nM final concentration. The hybridization buffer was composed of 20% Formamide, 2X SSC, 0.1g/ml Dextran Sulfate, 1mg/ml E.coli tRNA, and 2mM Vanadyl ribonucleoside complex, in nuclease free water. After overnight incubation in hybridization buffer and probes, cells were washed once in 20% Formamide and 2X SSC at 30°C for 30min, twice in 2X SSC at room temperature, stained with DAPI, and finally imaged in 2X SSC. smFISH imaging was performed on a Nikon Ti-E with Perfect Focus, Semrock FISH filtersets, Lumencor Sola illumination, 60x 1.4NA oil objective, and an Andor Zyla 4.2 sCMOS camera. Z-slices of DAPI, membrane-mTurquoise2 and smFISH were taken every 400nm through the sample. Segmentation of cellular boundaries was performed using the membrane targeted palmitoylated-mTurquoise2 with a 3D watershed algorithm. Dots were detected by thresholding on the distribution of local maxima of Laplacian-of-Gaussian kernel responses performed on each z-slice, with local-maxima defined around a 26-connected-pixel 3D region. Automatic image registration was performed in Matlab between the fluorescent protein in the final frame of the movie and DAPI stained image collected during smFISH imaging. RNA-seq On two separate days for biological replicates, ~500,000 were sorted of each subpopulation. Only the top 2% of reporters cells were collected in the positive gate for Zscan4c, while the lowest 50% were collected for the negative gate. For the Esrrb/Tbx3 double reporter (described above), 4% of the population made up the sorted double-negative population, 14% made up the sorted double positive population, and 63% made up the sorted Esrrb only population. The remaining unsorted cells made up buffer regions between subpopulations. Consistent with smFISH results, no Esrrb−/Tbx3+ population was observed. Differences in population fractions from smFISH-estimated population fractions is due to the half-life of the long-lived H2B-fused fluorescent proteins which only dilute by cellular division. Immediately following the sort, RNA was extracted using the Qiagen RNEasy Mini kit. 100 base single-end reads were generated on a HiSeq 2500. Galaxy was used to process RNAseq reads, using the Cufflinks package with default options. Briefly, reads were mapped using default Tophat parameters against the mouse mm10 genome. Cufflinks was used to estimate transcript abundance, and Cuffdiff was used to identify differentially expressed genes within E−T−, E+T−, and E+T+ sets, and then separately between Z+ and Z− sets.

QUANTIFICATION AND STATISTICAL ANALYSIS Kin Correlation Analysis

(KCA) applied to non-uniform cell state distributions In Box 1 of the main text, we derived a formula that related the two-cell correlation matrices to the transition matrix, under the assumption that the dynamics was reversible (or equivalently that it satisfied the condition of detailed balance) and briefly generalized the result to the case where all the states are not equally likely. Here, we will derive in full detail a formula for inferring the transition matrix from the observed two-cell correlation matrices assuming only that the dynamics is reversible. Then, in section 3 below, we will further relax the assumption of reversibility, and derive a more general expression using three-point correlation functions to infer transition dynamics. As in Box 1 of the main text, consider a transition matrix T ( I | M ) that represents the probability of observing a daughter cell in state I given that the parent cell was in state M . We assume that a given state I is observed in the population with frequency p I . With detailed-balance, for any given pair of states, the forward and reverse fluxes must be equal: T ( I ∣ M ) p M = T ( M ∣ I ) p I . The simpler condition used in Box 1 that T is a symmetric matrix, T ( I | J ) = T ( J | I ) is a special case of this expression, valid when p I = p J . A similar condition must hold going from a parent cell in state M to a descendent in state I after two generations: ∑ s T ( I ∣ S ) T ( S ∣ M ) p M = ∑ s T ( M ∣ S ) T ( S ∣ I ) p I , where s is summed over all possible state of the intermediate cell between the parent cell and its grand-daughter. The summation is equivalent to matrix multiplication and can be rewritten as, ∑ s T ( I ∣ S ) T ( S ∣ M ) = T 2 ( I ∣ M ) . More generally, for a cell in state M and its descendent u generations later in state I , the following condition must be satisfied: T u ( I ∣ M ) p M = T u ( M ∣ I ) p I . The joint probability of observing two cells at lineage distance u in states I and J is given by, C I J ( u ) = ∑ M T u ( I ∣ M ) T u ( J ∣ M ) p M . Reversibility of the dynamics implies that T u ( J | M ) p M = T u ( M | J ) p J . Making this substitution, we have, (1) C I J ( u ) = p J ∑ M T u ( I ∣ M ) T u ( M ∣ J ) = p J T 2 u ( I ∣ J ) To infer, we observe the correlation matrix C ( u ), and solve for the transition matrix T . Equation (1) can be rearranged to express the transition matrix in term of the correlation matrix, by first defining a rescaled correlation matrix, (2) C ∼ I J ( u ) = p J - 1 C I J ( u ) . It then follows that the transition matrix can be recovered by taking the appropriate root of the matrix C̃ , (3) T = C ∼ 1 / ( 2 u ) .

Kin Correlation Analysis

(KCA) applied to time-varying transition rates Previously, we assumed that transition rates remain constant over time. However, in a developmental context they could change systematically with time or generation number. In this subsection, we extend the above results to such cases. We still assume that the dynamics are stationary and reversible. As shown below, it is possible to fully recover time-varying dynamics by using the two-cell correlation functions at all lineage distances. Consider a time-varying transition matrix, T ( u ), where u denotes the number of generations back from the final time-point. u generations back, the probability of observing a daughter cell in state I conditional on the state M of its parent is given by the I , M th element of T ( u ). For an example of such dynamics, see Figure 5C in the main text. The two-cell correlation matrix for a pair of cells at lineage distance u takes the form, C I J ( u ) = p J ∑ M S u ( I ∣ M ) S u ( M ∣ J ) = p J S u 2 ( I ∣ J ) , where S is an effective transition matrix given by S u = T (1) T (2)··· T ( u ), and p J denotes the endpoint frequency of cells in state J . From a measurement of the two-cell correlation matrix C ( u ), S u can be inferred using Eqs. 2 and 3 , namely, define C ∼ I J ( u ) = p J - 1 C I J ( u ) . It follows, S u = C ∼ ( u ) To recover the time-varying transition rates T ( u ), we start at u = 1 and work our way backwards to larger values of u . The transition rates are given by, T ( 1 ) = S 1 , T ( 2 ) = T - 1 ( 1 ) S 2 , T ( 3 ) = T - 1 ( 2 ) T - 1 ( 1 ) S 3 , ⋮ T ( u ) = T - 1 ( u - 1 ) ⋯ T - 1 ( 2 ) T - 1 ( 1 ) S u . Lastly, we note that the above framework can be in principle extended to the case of continuous dynamics where the transition matrix is a continuous function of absolute time back to the common ancestor, i.e. where transitions have a probability per unit time (rather than per generation) of occurring, and where this probability itself changes with absolute time. To do so, the effective transition matrix, S u , is computed by taking the product integral of the continuous transition rate matrix, T́ , from the current time, t = 0, back to the time of the common ancestor, t = t c , namely, S ˋ u = exp ( ∫ t = 0 t = t c ln T ˋ ( t ) d t ) . This formulation accounts for biologically relevant cases in which the durations of cell cycles vary from cell to cell and/or over time.

Kin Correlation Analysis

(KCA) with probabilistic state assignment We show that the distribution of the Esrrb transcript counts in the low (E−) and high (E+) states overlapped significantly ( Fig. 2C ), such that it was not possible to assign a definite Esrrb state (either E− or E+) to a cell given a readout of its Esrrb transcript count. Here, we explain how we assigned probabilistic Esrrb states to each cell, and how the KCA framework is applied to probabilistic states. First, we fit a sum of two negative binomial distributions to the distribution of Esrrb transcript counts in single cells (black lines in Fig. 2C ). Let’s denote the distribution of transcript counts of the E− state as D − ( x ) and that of the E+ state as D + ( x ), where x is an integer denoting the transcript count in a given cell. More specifically, D − ( x ) is the probability that a cell in the E− state will have x Esrrb transcripts. The fit also has a free parameter that reflects the population fraction of each state. We will denote the population fraction of the E− state as f − and the population fraction of the E+ state as f + . It follows that f − + f + = 1. The probability that a cell with x Esrrb transcripts is in the E+ state is given by, p + = f + D + ( x ) f - D - ( x ) + f + D + ( x ) The probability that the cell is in the E− state is simply p − = 1 − p + . For large transcript counts, e.g. x = 200, f − D − ( x ) ≈ 0, which implies, p + ≈ 1. Alternatively, for some intermediate values of transcript counts, e.g. x = 75, f − D − ( x ) = f + D + ( x ), which implies p + ≈ 0.5, or that the cell is equally likely to be in the E− or the E+ state. Correlation matrices can be computed using probabilistic states in a similar manner as with definite states. However, whereas with definite state assignments, each pair of cells in states I and J contributes 1 to the I,J th element of the correlation matrix and 0 to all the other elements, with probabilistic state assignments, each pair of cells contributes P I P J to element C IJ of the correlation matrix. Lastly, the switching rates can be inferred from the correlation matrices computed using probabilistic states in a similar way as outlined in the previous section. However, in cases where the overlap between the two distributions is not symmetric, i.e. when it is more likely to misclassify a cell in the E− state as E+ than vice-versa, we need to first adjust the correlation matrices for incorrect assignment of states. On average, the probability that a cell in state J is assigned to state I is given by, Q I J = ∫ x f I D I ( x ) ∑ K f K D K ( x ) D J ( x ) d x where the integration runs over all possible values of transcript counts, x , and the summation K is over all states. Q is effectively a transition matrix satisfying the same properties as T ; for example, columns of Q sum to 1. However, unlike T , Q does not capture actual cell state transitions, but rather effective state transitions caused by measurement errors (for example, ambiguous mapping from transcript counts to cell state). Thus, we can imagine the dynamics as follows: the state of a cell u generations after its ancestor is given by the appropriate power of the transition matrix, namely, T u . The measurement error, at the endpoint, results in one additional mixing of states as given by matrix Q . Put together, the probability of observing a given state conditional on the state of the ancestor is given by the matrix QT u . To infer the actual transition matrix T , we must first remove the contribution of Q . The actual population fraction of the states, p̄ , can be calculated from the measured population fractions, p , as follows, p ¯ M = Q ^ M N p N where Q̂ is the inverse of the matrix Q . Similarly, the actual correlation matrix can be calculated from the measured correlation matrix as follows C ¯ = Q ^ C Q T ^ , where Q T ^ denotes the inverse of the transpose of matrix Q . The corrected populations fractions, p̄ , and correlation matrix, C̄ , can be used directly in Equations 2 and 3 above instead of p and C to infer the actual transition matrix. Computing the three-cell correlation functions for a general transition matrix with irreversible dynamics Here, we derive the general expression for the three-point correlation functions in terms of the transition matrix. As in the previous section, consider a transition matrix T ( I | M ) that represents the probability of observing a daughter cell in state I given that the parent cell was in state M . Unlike the previous section, we do not require that T ( I | M ) satisfies the condition of detailed balance, enabling analysis of cell state transition networks containing irreversible transitions. We would like to calculate the joint probability of observing three cells in states I , J , and K . The degree of relatedness of three cells is characterized by two lineage distance: u , the number of generations back to the common ancestor of the two more closely related pair of cells (observed to be in states J and K ), and v , the number of generations back to the common ancestor of all three cells (see Box 1 for a schematic). The three-cell correlation function takes the form, (4) C IJK ( u , v ) = ∑ M ( ∑ S T u ( K ∣ S ) T u ( J ∣ S ) T v - u ( S ∣ M ) ) T v ( I ∣ M ) p M , where the summation over S is over all possible states of the common ancestors of the two cells at lineage distance u . The summation over M is over all the possible states of the common ancestor of the three cells. p M is the expected probability of observing the common ancestor of all three cells in state M . For non-stationary dynamics, the probability of observing the common ancestor in a given state p M changes from generation to generation. However, p M is still related to the transition matrix in a self-consistent way. Namely, the probability that a cell u generations back will be in state M is given by (5) p M ( u ) = T u 0 - u ( M ∣ N ) p N ( u 0 ) where T u 0 − u ( M | N ) denotes element M , N of the transition matrix taken to the power of u 0 − u. u 0 is the number of generations back to the root of the tree, which is in state N with probability p N ( u 0 ). Equation (5) captures how the population fraction of each state changes over time as a function of the initial distribution of the states and the transition matrix. Equations (2) and (3) are general and do not require reversible (detailed balance) or stationary dynamics. Although an analytical solution for T in terms of C IJK ( u , v ) is not possible, we can solve the inference problem by considering the elements of the transition matrix as fitting parameters. We then calculate the expected three-cell correlation functions ( Eq. 4 ) and fit them to the observed three-cell correlation functions (see Methods in the main text for the numerical implementation). Simulating KCA for various types of dynamics Using KCA, we were able to accurately infer the underlying cell state transition network and transition rates in simulations by observing 30 cell pedigrees of 5 generations ( Figure S1 ). For reversible dynamics like those shown in Figure S1A,B , the transition network was inferred from the two-cell correlation functions. For networks with irreversible dynamics, like those shown in Figures S1C,D , we used the three-cell correlations for the inference. For example, in Figures S1B and S1C , which differ only in the reversibility of their transitions, the two-cell correlation functions are identical, but the three-cell correlation functions are different and can be used to infer the directionality of the transitions. Furthermore, we analyzed a previously published model of a 3-state system containing irreversible transitions among cancer cell states ( Gupta et al., 2011 ), and verified that KCA with three-cell correlations could indeed infer the previously determined rates ( Figure S1E ). Finally, we asked whether the KCA framework could be applied to non-stationary, branching cell fate determination networks, similar to those frequently observed in development and immunology. We simulated a 3-level branched cell fate tree with specific transition rates, applied KCA, and recovered the correct rates within statistical error ( Figure S1D ). This indicates that accurate inference is possible for branching fate trees with feasible amounts of experimental data. Here, we describe the details of the simulations used to generate the results shown in Figure S1 . Simulated pedigrees of 5 generations each were generated using Matlab. For S1A–C, the state of the root was selected randomly from the stationary distribution of cell states. In Figure S1D , the root was always set to the green state, resulting in a non-stationary distribution of cell states over the generations. At each generation, every node gave rise to two daughter nodes, whose states were selected randomly and independently from the probability distribution set by the state of the parent and the transition matrix. The two-cell and three-cell correlations were directly computed from the simulated pedigrees by measuring the frequency of occurrence of pairs and triplets of cell states at a given lineage distance over all pedigrees. We simulated 30 pedigrees for the plots in Figure S1A to C . For Figure S1D , we simulated 100 pedigrees. KCA using two-cell correlation functions was conducted on simulated data as outlined above using the framework in Box 1 and Supporting Information without any fitting. To infer the rates for irreversible dynamics, we used a set of fitting parameters, corresponding to the independent entries of a general asymmetric transition matrix. We then fit the three-cell correlation functions predicted from this transition matrix (see STAR Methods) to the observed three-cell correlation functions. A non-linear least square fitting algorithm was used (implemented in Matlab) to minimize the residual. Direct measurement of Esrrb switching dynamics We tracked and segmented each cell in time-lapse movies of colony growth using automated software and manual corrections, similar to previously described ( Singer et al., 2014 ). By integrating the background corrected pixel intensity in the nucleus of each cell, we obtained the accumulated level of H2B-mCitrine fluorescence at every point along each cell cycle. The rate at which fluorescence accumulated was used to estimate the promoter activity of Esrrb . To identify changes in promoter activity that corresponded to state switching, we fit either a single line, or two piece-wise linear segments to the fluorescence read-out of each cell using a least-squares method, implemented in Matlab. The first and last hour of each cell cycle was discarded to ensure reliable fluorescence read-out despite cell division. We used two criteria to identify state switching events: 1) the change in the slope across a division or between the segments of the two-line fit had to exceed a significance threshold. 2) A significant change in the slope (increase or decrease) had to persist into the subsequent cell cycle after division. The candidate switching events were identified automatically using a script implemented in Matlab and then verified manually. Computing the predicted spatial correlation functions for Esrrb We calculated the expected correlation of Esrrb state as a function of spatial separation distance from the inferred switching rates of Esrrb and the observed pedigrees as follows, C I J ( r ) = ∑ u q ( r ∣ u ) p ( u ) ∑ M T u ( I ∣ M ) T u ( J ∣ M ) p M where q ( r | u ) is the empirically determined probability of observing two cells at lineage distance u at spatial separation distance r; q ( r | u ) is plotted in Fig. 2Gii . p ( u ) is the probability that two randomly chosen cells will be at lineage distance u . This was empirically computed using the set of observed pedigrees. The expected and directly observed spatial correlation are plotted in Fig 2Gi . Selecting marker genes for the ES pluripotency states Previous studies have revealed that ES colonies exhibit a heterogeneous set of states potentially related to early embryonic cell types. For example, recent evidence identified a subpopulation of cells that express Zscan4, potentially corresponding to the totipotent 2 cell (2C)-state ( Falco et al., 2007 ; Macfarlan et al., 2012 ). This state is also associated with telomere-elongation, essential for long-term culture in vitro ( Zalzman et al., 2010 ) (although Zscan4 has also been shown to be activated by DNA damage responses and PI3K signaling (Storm et al., 2014)). Furthermore, representing slightly later stages of development, both inner cell mass (ICM)-like and epiblast-like stages can be identified and distinguished in culture by the high or low expression, respectively, of a cluster of correlated genes that includes Rex1, Nanog, and Esrrb . The totipotent state, marked by Zscan4 expression, shows low Rex1/Nanog/Esrrb expression ( Singer et al., 2014 ), potentially defining a sub-population among Rex1/Nanog/Esrrb -low cells. Finally, we identified a complementary sub-population within the Rex1/Nanog/Esrrb -high population, marked by expression of Tbx3. Tbx3 has been shown to destabilize pluripotency when lost or over-expressed. It also appears critical for mesendoderm specification, and its expression may change the global levels of DNA methylation in mouse ES cells (Dan et al., 2013; Ivanova et al., 2006 ; Lu et al., 2011 ; Niwa et al., 2009 ; Weidgang et al., 2013 ). However, it remains unclear how Tbx3 expression emerges dynamically from these states. As described in the main text, to verify that changes in the expression levels of the marker genes corresponded to collective changes in expression levels of multiple genes, we performed RNA-seq on subpopulations of cells sorted using fluorescent reporters for the three marker genes described above. We sorted out Esrrb/Tbx3 negative (E−T−), Esrrb-positive/Tbx3-negative (E+T−), and Esrrb/Tbx3 positive (E+T+) cells, as well as Zscan4 -positive (Z+) and –negative (Z−) cells, and observed hundreds of genes that were differentially expressed between these states ( Fig. 3D ), indicating that variations in marker gene expression do not simply reflect intrinsic noise, or fluctuations in the expression of individual genes, but rather indicate broad transcriptional changes. More specifically, compared with E+T− cells, E−T− cells expressed lower levels of pluripotency regulators ( Fig 3Di ), and higher levels of differentiation markers and signaling proteins ( Fig 3Dii ). In contrast, E+T+ cells showed reduced expression levels of signaling proteins and differentiation markers and increased levels of pluripotency genes compared to E+T−, suggesting that Tbx3 could mark a more pluripotent state. Moreover, we observed increased expression levels of 2C-associated genes like Tmem92, Tcstv3, Tdpoz3/4, and Zfp352 in the Zscan4 -positive cells compared with the Zscan4 -negative cells ( Fig 3Diii ). This result is consistent with Zscan4 marking the previously reported 2C-like state ( Macfarlan et al., 2012 ). Assigning cells to the pluripotent states The probabilistic assignment of the Esrrb state is presented in Figure 2C and the STAR Methods. T+ state was defined as Tbx3 transcript counts larger than 15. A threshold was obtained by comparing transcript counts to the direct observation of the promoter activity of Tbx3 gene in individual ES cells that had a knock-in fluorescent reporter for both endogenous loci of the Tbx3 gene (see Fig S2Ci, S2D ). Zscan4 expression levels were observed to be largely binary ( Figure 3A ). We used a threshold of 50 transcripts to assign cells to the Z+ state. Cells that were in the Z− and T+ states but were also in the E− state with a confidence level of at least 80% were assigned to the E−T+Z− state. We discarded any cells that were in the Z+ state but were also in the T+ state and/or the E+ state with a confidence level of 80% or higher (composing