The Protective Role of Symmetric Stem Cell Division on the Accumulation of Heritable Damage

Stem cell divisions are either asymmetric—in which one daughter cell remains a stem cell and one does not—or symmetric, in which both daughter cells adopt the same fate, either stem or non-stem. Recent studies show that in many tissues operating under homeostatic conditions stem cell division patterns are strongly biased toward the symmetric outcome, raising the question of whether symmetry confers some benefit. Here, we show that symmetry, via extinction of damaged stem-cell clones, reduces the lifetime risk of accumulating phenotypically silent heritable damage (mutations or aberrant epigenetic changes) in individual stem cells. This effect is greatest in rapidly cycling tissues subject to accelerating rates of damage accumulation over time, a scenario that describes the progression of many cancers. A decrease in the rate of cellular damage accumulation may be an important factor favoring symmetric patterns of stem cell division.

1 "Branching" model of ordered mutation accumulation We first consider the case in which mutations accumulate in a defined temporal order at a set of loci as shown in Fig. 1E and later (in Section 3) relax this assumption so that mutations may accumulate in any order.

Formulation
Let S represent a stem cell and D a non-stem cell. Assume the stem cell divisions are synchronized with a fixed cell-cycle duration of unity. When a stem cell divides, it spawns two daughter cells in one of three different ways with probabilities p r , p a and p e : Here, the fraction of divisions that are symmetric -producing daughter cells with a common fate -is denoted s (not to be confused with a selection coefficient). Numerical homeostasis of the stem-cell pool is ensured on average by balancing the probability of symmetric renewal, p r , with that of symmetric extinction, p e . Given a stage-i stem cell, S i , harboring neutral mutations at i = 0 . . . K − 1 loci, each of its daughters independently acquires a mutation at an additional locus, i + 1, with probability u i per stem-cell division. Thus, a symmetric renewal gives rise to one of four possible outcomes Similarly, an asymmetric division generates mutants according to whereas a symmetric differentiation looks like Since we restrict our exploration of parameter space to u i ≤ 10 −2 1, we may neglect probability contributions of order u 2 i (see Section 6 and Fig. S8). Making this approximation, and recognizing that the outcome X i + X i+1 (where X = S or D) cannot be distinguished from X i+1 + X i in a non-spatial model, we arrive at the following reduced set of outcomes Putting all this together, and book-keeping only the stem cells, we conclude that the division of a stem cell S i (i = 0 . . . K − 1) results in one of five possible outcomes with the following probabilities: This is a discrete-time multi-type branching process [1,2] in which each stem cell divides independently of all other stem cells. The process starts with N wild-type (stage-0) stem cells and ends when the first stage-K stem cell arises.

Efficient simulation algorithm
Let Z ij (t) be the (random) number of times that a category-j division in Eq. (S1) occurs at time t in a population of stage-i stem cells of (random) size N i (t). The stage sizes one cell cycle later are Notice that we have assumed that the stage-K stem cell has no birth-death dynamics. Instead S K → S K at each cell cycle. This simplification does not affect the conclusions we draw since the statistics that we compute relate to the appearance of the first stage-K stem cell and are thus independent of its dynamics. Conditioned on the sizes of the various stages, {N i (t) = n i }, the number of divisions in each category, Z i = (Z i,r , Z i,rm , Z i,a , Z i,am , Z i,e ), is distributed according to the multinomial distribution [3], where j z j = n i . This is an efficient way to simulate the stochastic process.

Deterministic analysis
Averaging the random variables in Eq. (S2) gives rise to terms of the form Z ij (t) that can be simplified by conditioning on the stage size as follows The factor in curly brackets is the average number of times that a category-j division occurs in a fixed-size population of stage-i stem cells and is given by the product of the stage size, n i , and the probability of a stem cell executing a category-j division, p ij . Making this replacement in the curly brackets immediately yields Using this result and the expressions for the category probabilities p ij given in Eq. (S1), one may then show that the average population sizes obey After a few cell cycles, t 1, these equations are well approximated by their continuous-time form If the mutation rates are identical and equal to u, then the relative abundance of the stages, N i /N , follows a truncated Poisson distribution We can also make analytical progress when the mutation rates are not all equal. Initially, the dynamics of each stage is determined solely by the influx from the previous stage, Using these solutions we may derive the following expression for the influx to stage i Plugging this expression into the influx-outflux inequality N i−1 u i−1 N i u i shows that the solutions in Eq. (S7) are valid provided which is met for most biologically realistic parameter values. Combining Eqs. (S8) and (S9) implies that the stage abundances are monotonically decreasing (S10)

Purely asymmetric risk
In a tissue undergoing a purely asymmetric pattern of division (s = 0), each stem cell behaves independently of all the others. This has two implications. First, the probability that the first K-fold mutant stem cell arises by time t, R K (t) = P [ T K < t ], is related to the corresponding risk per stem cell, r K (t), by where N is the population size. Second, the average number of stage-K stem cells is Putting these formulae together, we get (S11) Keeping N K constant while taking the limit N → ∞ yields This formula, together with Eqs. (S5), was used to screen for parameter sets with purely asymmetric lifetime risks in a defined range.

Heuristic analysis of protection
In this section we first derive the time scale on which a population of stem cells undergoing neutral drift is driven to extinction. We then use this result to heuristically derive the conditions under which symmetry is expected to delay mutation accumulation.

Median time for a stage to extinguish
Suppose that we observe n stem cells in a particular stage i of a purely symmetric population at time t = 0. Prior to mutation the stem-cell dynamics may be approximated by the continuous-time branching process, with average cell-cycle time of unity. The distribution of times, T e , at which the stage extinguishes is therefore [4,Eq. 8.58] where N i (t) is the random number of stage-i stem cells at time t. Scaling time according tô t = t/2n, the distribution takes on the form Keepingt constant while taking the limit n → ∞ yields the Fréchet distribution from Extreme Value Theory [5, p. 9], An extreme-value distribution is expected since T e is the maximum of the set of independent and identically distributed extinction times corresponding to the n lineages (clones) comprising the initial population. The Fréchet distribution has a rather heavy right tail, 1 − e −2n/t ∼ 2n/t as t → ∞, implying that the mean extinction time diverges logarithmically. We therefore used the median time to extinction, defined by P [ T e < τ e ] = 1 2 , and given by τ e = 2n ln 2 , (S14) as a measure of the time needed for a stage to extinguish given that it is observed to contain n stem cells. This expression is asymptotically exact at large stage sizes whereas it is a sufficiently good approximation for our purposes at small sizes. (The exact value implied by Eq. (S13) is 2 1−1/n / 1 − 2 −1/n .)

Protection criterion
We first determine which lineage (clone) of stage-i stem cells is responsible for progression to the next stage during the organism's lifetime, L. Eq. (S6b), which is valid provided L i/u i (see Eq. (S9)), may be re-cast as which says that the average lifetime abundance of stage-i stem cells approximately equals the average number of stage-i lineages founded during the course of life. In a symmetric population, most of the lineages never reach a substantial size and extinguish within a few cell cycles (Fig. 2D). Since each of these lineages has probability of order 1/x of reaching size x ≥ 1 or larger [6], on average only one of them will reach a size of at least N i (L) ≥ 1. Thus the mean lifetime abundance approximately measures the size of the largest clone likely to occur during a lifetime (asterisk in Fig. 2D). By virtue of its size, this clone carries the majority of the risk of progression to the next stage [6]. This intuitive picture fails when N i+1 (L) 1 because progression then occurs in rare lineages that grow to sizes much larger than N i (L) (see Section 2.3.1 and Fig. S2I).
Having identified the clone most susceptible to progression, we next determine the circumstances under which symmetric divisions delay its progression to the next stage. If the clone were dividing purely asymmetrically, it would progress in a time of order 1/ N i (L) u i . However, since it is dividing purely symmetrically, it will extinguish in a time of order N i (L) (derived in Section 1.5.1). Symmetric extinctions ought to delay progression provided the clone extinguishes faster than it progresses, N i (L) 1/ N i (L) u i . Thus the criterion for symmetry-dependent protection of stage i is and is valid provided L i/u i and N i+1 (L) ≥ 1 (notice that the last condition implies N i (L) ≥ 1 by virtue of Eq. (S10)). The threshold population size, 1/ √ u i , represents the size a symmetric stage-i lineage needs to reach in order to progress to the next stage [6]. Re-casting Eq. (S8) in the form and plugging this expression for the abundance into Eq. (S15), we arrive at an alternate form of the protection condition again valid for L i/u i and N i+1 (L) ≥ 1. Eqs. (S15) and (S17) suggest natural scales on which to measure abundance and mutation rate, respectively, Eq. (S16) may then be used to establish the following simple relation between the scaled parameters valid for L i/u i .

Probability that a lineage eventually mutates
Given that there is initially just one wild-type (stage-0) stem cell, we wish to calculate the probability that eventually one of its descendants mutates, where P 0 (t) is the probability of surviving the mutation by time t. Since the stage-1 stem cell has no dynamics (K = 1), the event (N 1 (t) = 0) is equivalent to the nonoccurrence of the mutation by time t, where N 0 (t) and N 1 (t) denote the numbers of wild-type and mutant stem cells, respectively. To carry out the calculation we first need to introduce the theory of multi-type branching processes. Using the shorthand (n 0 , n 1 ; t) to represent the event (N 0 (t) = n 0 , N 1 (t) = n 1 ), we define The probability generating functions, F 0 and F 1 , satisfy the recursion relations [ where we have dropped the arguments (x 0 , x 1 ) for clarity, and are the probability generating functions evaluated at the first division.
Having introduced the theory of multi-type branching processes, we now see that the probability that the lineage descending from a stage-0 stem cell survives a mutation is just P 0 (t) = F 0 (1, 0; t). Using the recursion relations, Eq. (S21), we deduce that where the second argument of f 0 is since the mutant cannot extinguish. We may rewrite the kinetic equation for P 0 (t) as (S22) The long-time limit, P (∞) 0 = lim t→∞ P 0 (t), therefore satisfies the quadratic equation The physically meaningful solution is given by the smallest root and can be written in the form P It is convenient in the analysis that follows to re-write the expression under the square root as Biologically plausible values for the mutation rate satisfy u 0 1, in which case β 0 ≈ 1 and If the symmetry fraction is small compared with the mutation rate, s u 0 , then γ 0 ≈ u 0 s , and the square root is approximated by 1 + u 0 s , so that the lineage mutation probability is In this regime, the lineage behaves as if it were dividing purely asymmetrically. Now consider the alternative limit where the symmetry fraction is large compared with the mutation rate, s u 0 . There are two sub-cases: s ∼ 1 and s 1. In both cases γ 0 is of the same order as 2u 0 /s so that the leading-order term in the lineage mutation probability is Generally, the probability that the lineage descending from a stage-i stem cell accumulates one further mutation is ( Fig. S2 Conceptually similar manipulations can be employed to calculate the probability that a lineage accumulates two or more mutations (see also [6]).

Mean time that a lineage drifts before mutating
Given that the lineage descending from a wild-type (stage-0) stem cell picks up a mutation, what is the mean time it drifts before doing so, τ 0 ? Remembering that stage-1 stem cells have no dynamics, we may write where T 1 is a random variable representing the time at which the first single-mutant stem cell arises in the lineage. The mean drift time τ 0 is related to the conditional probability of surviving the mutation, where the conditioning on N 0 (0) = 1, N 1 (0) = 0 has been made implicit in the interests of clarity. The summand is related to the unconditioned survival probability in Eq. (S20) by Using this relation in Eq. (S25), and replacing the sum by an integral (a valid approximation for large drift times, τ 0 1), the mean time until the lineage mutates becomes The survival probability P 0 (t) satisfies a Riccati equation, which is the continuous-time approximation to Eq. (S22). The time-dependent solution is and r are the two roots of Eq. (S23). As shown in Section 1.6, the roots may be approximated by where, as in Section 1.6, If the symmetry fraction is small, s u 0 , then P (∞) 0 ≈ 0 (see Section 1.6), r ≈ 2u 0 /s 1, and the survival probability is approximated by its purely asymmetric form, Performing the integral in Eq. (S26) yields which is indeed large, τ 0 1. If, on the other hand, the symmetry fraction is large, s u 0 , then (see Section 1.6). After some algebra, Eq. (S26) reduces to Performing the integral yields which is again large, τ 0 1. We note that at small times, t τ 0 , we have θ ≈ 1 − √ 2u 0 s t, so that the survival probability, Eq. (S27), reduces to P 0 (t) ≈ 1 − u 0 t, which, Eq. (S28) tells us, is also the short-time form in the limit of small symmetry fractions, s u 0 . In other words, the probability that a lineage mutates by time t τ 0 is given by irrespective of the value of the symmetry fraction. Generally, given that the lineage descending from a stage-i stem cell picks up an additional mutation, the mean time it takes to do so is ( Fig. S2(B))

"Moran" model of ordered mutation accumulation
The total population size fluctuates under the branching model. These fluctuations are relatively insignificant in large populations but can significantly violate numerical homeostasis in small populations. Yet real tissues attenuate errant fluctuations via feedback mechanisms. We therefore developed a phenomenological model of feedback that captures all of the stem cell divisions defined by Eq. (S1). This model is a variation on the well-known Moran model from population genetics [9; 10, p. 79] in which time is continuous and the cell cycle period is exponentially distributed with a mean value of unity.

Formulation
We define first how the asymmetric divisions are treated in the model. The rate at which asymmetric divisions occur in a population of stem cells is N (1 − s), where N is the timeinvariant total number of stem cells and s is the probability that a stem cell division is symmetric. When such a division occurs, it occurs in a stage-i stem cell with probability n i /N , where n i is the number of stage-i stem cells. Only asymmetric divisions in which the daughter mutates need to be simulated since only those change the state of the populationdecrementing the number of stage-i stem cells while incrementing the number of stage-(i + 1) stem cells. The net rates of such events are We turn next to the symmetric divisions. In the Moran approach, each populationincrementing symmetric renewal is accompanied by a population-decrementing symmetric extinction. Together these divisions constitute one replacement event. Balancing renewals with extinctions in this way ensures that the net population size is strictly constant. Since each replacement constitutes two stem cell divisions, replacements in a purely symmetric population must occur at half the rate of stem cell divisions in an equivalent purely asymmetric population. Thus the rate of replacements in a population undergoing a mix of replacements and asymmetric divisions is N 1 2 s. When a replacement occurs, a stage-i stem cell is chosen for renewal with probability n i /N whereas a stage-j stem cell is chosen for extinction with probability n j /N . The symmetric renewal results in mutation of one of the two stem cell daughters with probability 2u i , in which case the number of stage-(i + 1) stem cells is incremented whereas the number of stage-j stem cells is decremented. The case i + 1 = j need not be simulated since it does not change the state of the population. The net rates of such events are Alternatively, with probability 1 − 2u i , the symmetric renewal is of the form S i → S i + S i , in which case stage-i is incremented whereas stage-j is decremented. Again the case i = j need not be simulated. The net rates of these events are This process is simulated by choosing the time interval until the next event by sampling from an exponential distribution with intensity equal to and executing one of the events described above with probability proportional to its rate (the Gillespie algorithm). As expected, it coincides with the branching model at long times and large population sizes (Fig. S1A).

Deterministic analysis
We first derive the deterministic equations for a general continuous-time stochastic process. Letting ( n, t) denote the event ( N (t) = n), the expected size of the stage-i stem cell population is with transition probability p n→ n = P ( n , t + δt | n, t).
It is useful to decompose the term in curly brackets as follows Plugging this into Eq. (S34), dividing across by δt, and taking the limit δt → 0, we arrive at an equation for the mean size of the stage-i stem cell population where Having derived a general equation governing the time dependence of the mean population sizes, we next apply it to the stochastic process described in Section 2.1. Given a vector n = (n 0 , n 1 . . . n K ) describing the numbers of stage-0, 1 . . . K stem cells in a population, Table 1 classifies the states n = (n 0 , n 1 . . . n K ) that can be reached in a single transition and that have n i = n i . Inserting these state transitions into Eq. (S37) reduces it to The summation in Eq. (S39a) can be performed immediately Similarly, the summation in Eq. (S39b) reduces to σ (s) The summation in Eq. (S39c) is more involved but can be performed by partitioning its summands into three terms as follows The summation over a i,j can be performed immediately since n K = 0 during the course of the simulation. The summations over the terms involving b i,j can be written where the square bracket indicates that the corresponding term should be omitted from the sum. Most of these terms cancel when Eq. (S44) is subtracted from Eq. (S43) yielding Adding Eqs. (S42) and (S45) yields Plugging the expressions in Eqs. (S40), (S41) and (S46) into Eq. (S38) and performing some algebra yields Finally, plugging this expression into Eq. (S36) yields the deterministic equations for the feedback model These equations are identical to those derived in the branching model at long times, Eq. (S5), as they should be.

Formal analysis of the accumulation of two mutations
Previous studies in the field of population genetics have identified a number of parameter regimes in which mutation accumulation in populations subject to genetic drift can be described using simple analytical formulae [6,[11][12][13][14][15]. Here, we build on this theoretical work to better understand why symmetric stem cell divisions protect against mutation accumulation and to derive mathematical formulae that accurately describe the extent of protection.

Derivation of PF formula
The generation of the first double-mutant stem cell in the purely asymmetric trajectory of Fig. S2C can be described by a pair of independent exponential random variables, X 0 and X 1 , with rate parameters N u 0 and u 1 , respectively. These values for the rates follow from the fact that a population containing n i stage-i stem cells generates a stage-(i + 1) stem cell for the first time at rate where the rates λ (am) i and λ (sm) ij are defined in Section 2.1. The time at which the first double-mutant stem cell arises, X 0 + X 1 , is hypoexponentially distributed [16, p253] with cumulative distribution R 2 given by and initial condition R 0 (0) = 1, R 1 (0) = 0 and R 2 (0) = 0. The solution is (Fig. 3E) where the 2-parameter hypoexponential distribution is given by This expression for R 2 is a good approximation to the exact value provided the mean abundance of single-mutant stem cells is small, N 1 (t) 1, a regime we call the "Stochastic" regime. In the alternative "Deterministic" regime, N 1 (t) 1, the single-mutant abundance is well approximated by its mean, N 1 (t) ≈ N 1 (t) , and the cumulative risk of generating a double-mutant stem cell is which is a special case of Eq. (S12). To progress further analytically, we solve Eq. (S47) to get an expression for the mean single-mutant abundance which can be approximated by provided there is a separation of timescales, u 0 u 1 . We needn't consider further the regime t 1/u 0 because we carried out our simulations at u 0 = 10 −6 whereas biologically relevant timescales are no larger than ∼ 10 4 cell cycles. The boundary between the Stochastic and Deterministic regimes, N 1 (t) ∼ 1, is therefore defined by (Fig. 3I) Plugging Eq. (S52) into Eq. (S50), one may show that the cumulative risk in the Deterministic regime, N u 1 /u 0 , t 1/N u 0 , is Let us now consider the purely symmetric trajectory in Fig. S2D. Many single-mutant lineages are expelled from the tissue before one arises that survives neutral drift long enough to acquire a second mutation. This occurs with probability Q (∞) 1 = √ 2u 1 (see Eq. (S24)), and when it does, the lineage drifts for a time τ 1 = 2 ln 2/ √ 2u 1 before mutating (see Eq. (S30); inset of Fig. S2D). The single-mutant lineage destined to mutate typically does not fix in the population -an adaptation mechanism known as "Stochastic Tunneling". In this regime, the purely symmetric cumulative risk is (Fig. 3E) fails to protect when the tissue age is much smaller than the drift time, t τ 1 (Fig. 3E). Here, the second mutation typically occurs in the first single-mutant lineage (Fig. S2E, F) because, for this particular parameter set, single-mutants rarely arise so early, Thus T 1 , the time at which the first single-mutant stem cell arose, has probability density function Moreover Eq. (S29) tells us that, if the single-mutant lineage is founded at time t , then it yields a double mutant by time t with probability where T 2 is the time at which the first double-mutant stem cell arose. The net cumulative risk of acquiring the double-mutant stem cell by time t is therefore (Fig. 3E) irrespective of the symmetry fraction. This expression coincides with the purely asymmetric risk derived in Eq. (S48) since t min{1/N u 0 , 1/u 1 } (for the parameters of Fig. S2C-F). In this "short-time" regime (t τ 1 , Fig. 3I), the beneficial effect of symmetric extinctions are cancelled by the deleterious expansion of single-mutant clones (Fig. 1C, D). Fig. S2G-I shows that the short-time regime can extend to tissue ages as large as 10 3 stem cell cycles if the secondary mutation rate is small enough, u 1 1/t 2 , explaining why protection vanishes in the lower part of Fig. 3D. Stochastic tunneling occurs for population sizes up to N = 1/u 0 , beyond which the symmetric single-mutant abundance is approximately deterministic [6,13] and protection is lost (see Section 2.3.2).
When the population size is small, the probability that a single-mutant clone fixes in the population exceeds the probability that it mutates before fixing, 1/N Q (∞) 1 (Fig. 3I). This is the "Sequential Fixation" regime [6,13] where it becomes possible that the first doublemutant stem cell arises after fixation of the corresponding single-mutant lineage. To make analytic progress in this regime, we consider only sufficiently large times, t τ 1 , allowing us to neglect the fixation time (which is of order N τ 1 ). The probability that a single-mutant lineage founded at time t = 0 mutates by some later time t via sequential fixation is whereas via stochastic tunneling it is Using the fact that τ 1 1/N u 1 , one may then show that the sequential fixation route is more likely than the stochastic tunneling route, P [ T 2 < t , fix ] P [ T 2 < t , no fix ]. At long times (Fig. 3I), one may neglect the time between (instantaneous) fixation of the single-mutant lineage and acquisition of the second mutation (Fig. S2L). Thus the cumulative symmetric risk in this regime is determined by the rate at which single mutants destined for fixation arise. This rate is u 0 (single-mutant lineages arise at rate N u 0 and fix with probability 1/N ) implying that the net cumulative risk is (Fig. S2J) On long time scales, t 1/N u 1 , symmetric stem cell extinctions flush many single-mutant lineages out of the tissue until one survives drift, fixes, and mutates leading to significant protection (Fig. S2M). At short times, τ 1 t 1/N u 1 ≤ 1/N u 0 (we assume that u 1 ≥ u 0 ), the second mutation typically occurs in the first single-mutant lineage (Fig. S2O). Furthermore, Eq. (S57) implies that We may therefore proceed exactly as we did for the Stochastic Tunneling regime at short times obtaining Eq. (S56), which again coincides with the purely asymmetric risk derived in Eq. (S48) since t 1/N u 1 ≤ min{1/N u 0 , 1/u 1 }. The beneficial effect of symmetric extinctions is thus cancelled by the deleterious effect of clonal expansion, as we found in the Stochastic Tunneling regime at short times.
In summary, symmetry is protective provided the tissue cycles rapidly and/or the secondary mutation rate is high (Fig. S3A-C shows that symmetry is not necessarily protective if instead the primary mutation rate is fast). Together, Eqs. (S48), (S54), (S55) and (S59) can be used to estimate PF = R A /R S to an accuracy of 40% (Fig. S2P-R) throughout the protected part of parameter space (Figs. 3I-L). Fig. S3B shows that when the average number of new single mutants produced per cell cycle becomes large, N u 0 1, the purely symmetric single-mutant abundance is well approximated by its mean. This regime is included in the Deterministic regime of a purely asymmetric population, defined by Eq. (S53) (see also Fig. 3I), and so the cumulative risks coincide (Fig. S3C). When N u 0 1, the time at which the first successful single-mutant lineage appears is small compared to its drift time, 1/N u 0 Q (∞) 1 τ 1 . In other words, the rate of progression is limited by the time taken for a successful single-mutant lineage to drift before mutating. This assumes however that the first such lineage gives rise to the first double-mutant stem cell. The fact that the mean first passage time to the double mutant (approximately the time at which cumulative risk reaches 50% in Fig. S3C) is significantly smaller than τ 1 indicates instead that many single-mutant lineages founded early on compete against one another until one acquires the first double-mutant stem cell improbably early. We conclude that, in the Deterministic regime of a symmetric population, protective clonal extinctions are out-striped by the rapid production of single-mutant stem cells.

Unordered accumulation of two mutations
In this section we generalize the models of mutation accumulation presented in Sections 1 and 2 to include the possibility that loci mutate in any order. We consider only the case of two loci, A and B.

"Branching" model
Assume stem cell divisions are synchronized. A wild-type cell, X 0 (where X is either a stem cell S or non-stem cell D), acquires a mutation at locus A with probability u 0→aB and at locus B with probability u 0→Ab per stem-cell division, independent of the fate of its sister cell. Each of these single mutants, X i∈{aB,Ab} , becomes a double-mutant, X 2 , with probability u i→2 per stem cell division. All other transitions (e.g. 0 → 2) are not permitted in a single stem-cell division (u 0→2 = 0). Following Section 1.1, one may show that the division of a wild-type stem cell results in one of seven possible outcomes with the following probabilities: Similarly, the division of a single-mutant stem cell with genotype i ∈ {aB, Ab} results in one of five possible outcomes: The update rules for the random population sizes, N i , of each of the four possible genotypes are i ∈ {aB, Ab} and Z i∈{0,aB,Ab} are drawn from the multinomial distribution, Eq. (S3).

"Moran" model
We formulate the model analogously to Section 2.1. The net rate of asymmetric divisions that decrement the number of stem cells with genotype i ∈ {0, aB, Ab}, n i , while incrementing the number of stem cells with genotype k ∈ {aB, Ab, 2} is where s is the probability that a stem cell division is symmetric. There are two types of replacements. In the first type, which occurs at rate a stem cell with genotype i ∈ {0, aB, Ab} is chosen for symmetric renewal, a stem cell with genotype j ∈ {0, aB, Ab} is chosen for symmetric extinction, and a stem cell with genotype k ∈ {aB, Ab, 2} is produced by mutation of one of the daughters of the symmetric renewal. This type of replacement thus increments the k-population and decrements the j-population such that the total number of stem cells, N , remains unchanged. Alternatively, with rate λ (s) the symmetric renewal is of the form S i → S i + S i , in which case the i-population is incremented whereas the j-population is decremented.

Deterministic solution
Following Section 1.3, one may show that the average population sizes obey The solution of these equations is where H 2 is the 2-parameter hypoexponential distribution defined in Eq. (S49).

Purely asymmetric risk
Eq. (S11) tells us that the probability that the first double-mutant stem cell arises by time t is

Purely symmetric risk when loci mutate independently
Consider the case depicted in Fig. 5B where t τ aB = 2 ln 2 √ 2u aB→2 s . Some double mutants arise from aB single mutants via stochastic tunneling whereas others arise from a deterministic background of Ab single mutants. The net rate of mutation accumulation is approximately the sum of two independent Poisson processes with rates 2u aB→2 , and k Ab (t) = N Ab (t) u Ab→2 , respectively [6]. For the parameter values of Fig. 5B, and the cumulative risk of generating a double-mutant stem cell is To combine mutational paths at short times, t τ aB , it is helpful to introduce a random variable Y 1 ∈ {S aB , S Ab } representing the genotype of the single-mutant stem cell that generates the first double-mutant stem cell. In cases where an aB stem cell generates the first double mutant (Y 1 = S aB ), typically it is the first aB stem cell that does so (Fig. S3E). Therefore, we may invoke the argument leading to Eq. (S56) obtaining In the alternate scenario -the double mutant arises from a deterministic background of Ab single mutants -the cumulative risk is for the parameter values of Fig. S3D-H. The mutational paths are combined by simply adding the cumulative risks in Eqs. (S69) and (S70) obtaining ( Fig. S3H)

Effect of selection on ordered mutation accumulation
In our analysis thus far, we assumed each stem cell contributes on average one stem cell to the population one cell cycle later. This is reasonable for asymmetric stem cell divisions but for symmetric divisions, renewals and symmetric extinctions could be imbalanced (ie. subject to selection). Here, we relax the assumption of neutrality by allowing the mean number of stem cell descendants contributed by a given parent to depend on its stage, subject to the constraint that the net stem cell population size is conserved.

Formulation
Leaving the asymmetric divisions in the Moran model of Section 2.1 unchanged (i.e. retaining the rates in Eq. (S31)), we generalize the model such that, when a stem cell replacement occurs, a stage-i stem cell is chosen for symmetric renewal with probability w i k w k n k , where w i is the fitness of a stage-i stem cell and n i is the number of such stem cells. Thus, the rates of replacements become where i, j = 0 . . . K − 1 index the stages, s is the fraction of stem cell divisions that are symmetric, and is the mean fitness of the population. Notice that when all stages have equal fitness the transition rates in Eq. (S72) reduce to Eqs. (S32) and (S33) for a neutral stem cell population.

Deterministic Analysis
In this section, we generalize the deterministic equations, Eq. (S47), to incorporate selection. As shown in Section 2.2, the mean number of stage-i stem cells, N i , evolves according to where n = (n 0 . . . n K ) represents the numbers of stage-0 . . . K stem cells in a population and P ( n, t) is the probability of finding the population in the configuration n at time t.
Repeating the derivation of Λ (i) n presented in Section 2.2, this time using the replacement rates defined by Eq. (S72), we obtain In contrast to the neutral case, Eqs. (S73) are no longer closed and depend upon factors of the form N i /w , for which one must derive further differential equations. We can truncate this hierarchy by making the approximation N ≈ N leading to where the "ensemble average" of the mean population fitness is This "mean-field approximation" is expected to fail whenever one or more stages are small, under which circumstances the distribution P ( n, t) becomes very broad and is poorly approximated by its mean.

Protection in the compartmentalized intestine
In an intestine undergoing a purely asymmetric pattern of division, each stem cell behaves independently of all the others. Thus, as argued in Section 1.4, the intestine-wide risk of ordered accumulation of K mutations with mutation rates u 0 , . . . , u K−1 is (red lines in Fig. S7) where N is the number of stem cells per crypt and M is the number of crypts. The function H K (k 0 , . . . , k K−1 ; t) is the K-parameter hypoexponential distribution [16, p253], which we computed by solving the ordinary differential equations with initial conditions P 0 = 1, P i = 0 (i = 1, . . . , K − 1), H K = 0. In the purely symmetric case (under the Moran model; Section 2.1), the K-fold mutant arises in an individual crypt via sequential fixations (provided the number of stem cells per crypt and the mutation rates are small enough [6,17]) implying that the intestine-wide risk is (green lines in Fig. S7) 6 Effect of allowing for simultaneous mutations in both daughter stem cells In Section 1.1, we formulated the "Branching" model of ordered mutation accumulation by neglecting the possibility that, upon division of a stem cell, both daughter cells may simultaneously mutate. To find out whether this approximation is good, we re-formulate the model to include this possibility. Retaining probability contributions of order u 2 i , one finds that the division of a stem cell S i (i = 0 . . . K − 1) now results in one of six possible outcomes with the following probabilities: Let Z ij (t) be the (random) number of times that a category-j division in Eq. (S77) occurs at time t in a population of stage-i stem cells of (random) size N i (t). The update rules for the N i are again given by Eq. (S2), where this time but Z i is again sampled from the multinomial distribution in Eq. (S3). The recursive equations for the average population sizes again follow Eq. (S4). We re-ran stochastic simulations using the more exact algorithm formulated here; the results are unchanged (compare Fig. S8 with Fig. 1G).   Table S1 were classified into 4 types based on the number of stochastic stages. Representative symmetric trajectories are shown. Notice the correlation between the number of stochastic stages and mean PF (averaged over all parameter sets with a given number of stochastic stages). 1/N u 1 , singlemutant lineages frequently extinguish before one survives drift, fixes in the population, and then rapidly acquires the next mutation. The time taken to fix, N , and to acquire the second mutation once fixed, 1/N u 1 , are negligible compared to the time taken for a singlemutant lineage destined for fixation to arise, 1/u 0 . (O) At short times, t 1/N u 1 , the time between fixation and mutation cannot be neglected. (M) Symmetric extinctions outcompete fixation events to reduce mutation accumulation risk. (P -R) Accuracy of piecewise analytic formula for PF measured by the ratio of analytically computed PF (Fig. 3J-L) to exact PF as calculated via Monte Carlo simulation (Fig. 3F-H). Grey contours delineate regions (red) where the fractional error of the analytic formulae |PF analytic − PF exact | /PF exact is less than 40%, including practically all the protected zone (PF exact > 2; white contour). Panels A and B were generated under the Branching model, Eq. (S1), whereas panels C -R were generated using the Moran model, Section 2.1. In panels C -O, purely asymmetric (symmetric) trajectories are in red (blue). cc, cell cycles. Figure S3: Unordered "fast" and "slow" mutations (A -C) A "fast-slow" ordered pathway. (A, B) The abundance of Ab stem cells is approximated by its mean value, N 1 (t) ≈ N (1 − e −u 0 t ), which follows from Eq. (S51) when t, 1/u 0 1/u 1 . (C) Simulated cumulative risk (symbols) is approximated by 1 − exp − t 0 N 1 (t ) u 1 dt . (D -H) Dynamics at short times, t τ aB , of unordered "fast" and "slow" loci. The black line in panels F and G is Eq. (S67). (H) Simulated cumulative risk (symbols) is approximated by Eq. (S71). In all panels, population size is N = 10 4 stem cells and mutation rates are 10 −2 ("fast") and 10 −6 ("slow"). Time courses are plotted until the first double-mutant stem cell appears in the entire stem cell population. The purely asymmetric trajectory (s = 0%) is representative of all trajectories examined whereas the mixed (s=10%) and purely symmetric (s=100%) trajectories show the most frequently observed type since all four possible combinations of stochastic tunneling and sequential fixation were observed at appreciable frequencies (see also [17]). The intestine was assumed to comprise 10 6 crypts, each containing 10 stem cells. Mutation rates are 10 −6 (slow) and 10 −3 (fast).

Figure S5: Protection persists when selection acts on stochastic stages
Protection against ordered accumulation of K = 2 mutations after 1000 stem cell cycles for symmetry fractions s = 100% (A) and 10% (B), calculated by Monte Carlo simulation of the generalized model presented in Section 4.1. The selection coefficient is defined in the model by (w 1 − w 0 )/w 0 , where w i is the fitness of stage i (see Section 4.1). The insensitivity of PF to wide variations in the selection coefficient is an example of the general principle in population genetics that selection is ineffective provided the magnitude of the selection coefficient is smaller than the inverse population size. Mutation rates are u 0 = 10 −6 and u 1 = 10 −3 per locus per stem cell cycle. Figure S6: "Increasing" mutation rates yield a broad distribution of latencies Probability distributions of times at which the first double-mutant stem cell arose in a population of N = 10 4 symmetrically dividing stem cells, from simulations of the discretetime branching process defined by Eq. (S1) (bars) and from the probability mass function P [ T 2 = t ] = k(t) exp − t 0 k(t ) dt (lines). In the Deterministic regime (A), the rate constant is given by k = N 1 u 1 , where the mean abundance of single-mutant stem cells is N 1 ≈ N (1 − e −u 0 t ), whereas it is k = N u 0 √ 2u 1 in the Stochastic Tunneling regime (B). When the mutation rates are decreasing, the distribution of latency until the first doublemutant stem cell is narrow (A), but when the mutation rates are increasing, the distribution becomes wider (B), even at the same mean latency (858 cell cycles in both cases). Histogram bar height represents the probability that the mutation occurred between the bar edges. Insets show typical stochastic realizations (red) and mean single-mutant abundance (black). Figure S7: Protection in the human colon Cumulative risk of ordered accumulation of K = 4 mutations by ages 30 (A, C) and 60 (B, D), assuming mutation rates increase 1000-fold (A, B) or 100-fold (C, D) during the course of mutation accumulation from an initial rate of u 0 = 5 × 10 −7 per locus per stem cell cycle [18]. Lines are Eqs. (S75) and (S76) whereas symbols are Monte Carlo simulations (under the Moran model; Section 2.1). The colon is assumed to be compartmentalized into M = 10 7 crypts, (∼10 4 crypts/cm 2 × ∼10 3 cm 2 /colon [19]) each containing N = 20 stem cells [20,21] dividing purely asymmetrically (red) or symmetrically (green) 100 times per year. Mutation rates of consecutive stages (u 0 , u 1 , u 2 , u 3 ) are (A) 5 × 10 −7 , 5 × 10 −7 , 5 × 10 −6 , 5 × 10 −4 ; (B) 5 × 10 −7 , 5 × 10 −7 , 5 × 10 −7 , 5 × 10 −4 ; (C) 5 × 10 −7 , 5 × 10 −7 , 5 × 10 −5 , 5 × 10 −5 ; (D) 5 × 10 −7 , 5 × 10 −7 , 5 × 10 −6 , 5 × 10 −5 . Figure S8: Symmetry protects even when mutations may occur simultaneously in both daughter stem cells The more accurate ordered mutation accumulation model presented in Section 6 was used to generate a distribution of PFs over a random ensemble of parameter sets equivalent to that used in Fig. 1G (Materials and Methods; Table S1). The distribution is unchanged (within sampling error).