
==== Front
J Math Biol
J Math Biol
Journal of Mathematical Biology
0303-6812
1432-1416
Springer Berlin Heidelberg Berlin/Heidelberg

39222150
2134
10.1007/s00285-024-02134-4
Article
Error-induced extinction in a multi-type critical birth–death process
http://orcid.org/0009-0004-6180-1685
Guasch Meritxell Brunet xell.brunetguasch@ed.ac.uk

1
Krapivsky P. L. 23
Antal Tibor 1
1 grid.4305.2 0000 0004 1936 7988 School of Mathematics and Maxwell Institute for Mathematical Sciences, University of Edinburgh, Edinburgh, EH9 3FD UK
2 https://ror.org/05qwgg493 grid.189504.1 0000 0004 1936 7558 Department of Physics, Boston University, Boston, MA 02215 USA
3 https://ror.org/01arysc35 grid.209665.e 0000 0001 1941 1940 Santa Fe Institute, Santa Fe, NM 87501 USA
2 9 2024
2 9 2024
2024
89 4 3626 11 2023
2 7 2024
9 8 2024
© The Author(s) 2024
2024
https://creativecommons.org/licenses/by/4.0/ Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Extreme mutation rates in microbes and cancer cells can result in error-induced extinction (EEX), where every descendant cell eventually acquires a lethal mutation. In this work, we investigate critical birth–death processes with n distinct types as a birth–death model of EEX in a growing population. Each type-i cell divides independently (i)→(i)+(i) or mutates (i)→(i+1) at the same rate. The total number of cells grows exponentially as a Yule process until a cell of type-n appears, which cell type can only divide or die at rate one. This makes the whole process critical and hence after the exponentially growing phase eventually all cells die with probability one. We present large-time asymptotic results for the general n-type critical birth–death process. We find that the mass function of the number of cells of type-k has algebraic and stationary tail (size)-1-χk, with χk=21-k, for k=2,⋯,n, in sharp contrast to the exponential tail of the first type. The same exponents describe the tail of the asymptotic survival probability (time)-ξk. We present applications of the results for studying extinction due to intolerable mutation rates in biological populations.

Centre for Doctoral Training in Mathematical Modelling, Analysis and Computation, University of Edinburgh (GB)EP/S023291/1 Guasch Meritxell Brunet issue-copyright-statement© Springer-Verlag GmbH Germany, part of Springer Nature 2024
==== Body
pmcIntroduction

Genetic defects in DNA replication fidelity and repair result in mutator phenotypes that accelerate the adaptation of microbes and cancer cells. However, certain combinations of mutator alleles increase mutation rates to intolerable levels, resulting in error-induced extinction (EEX), where every cell eventually acquires a lethal mutation (Morrison et al. 1993; Fijalkowska and Schaaper 1996). Error-induced extinction occurs within a few generations in bacteria and haploid yeast with mutations on both DNA proofreading and repair (Morrison et al. 1993; Fijalkowska and Schaaper 1996; Herr et al. 2011). Tumours with mutator phenotypes reach an upper limit of mutations, which has been interpreted as evidence for a maximal mutation rate in cancer (Fox and Loeb 2010; Schumacher et al. 2019).

Defining the maximum mutation rate is important for understanding the long-term fitness of mutator alleles, with applications for studying the evolution of hyper-mutated tumours and the synthetic lethality of mutator alleles (Topatana et al. 2020). Experimentally, this amounts to producing cell lines with mutator mutations, and measuring the mutation rates when colonies are viable. This strategy has identified intolerable mutation rates for bacteria (Morrison et al. 1993), yeast (Fijalkowska and Schaaper 1996; Soriano et al. 2021) and cancer in mice (Albertson et al. 2009). However, beyond a threshold, the cell lines are not viable or acquire antimutator alleles. Thus, it is not possible to study the behaviour of populations at the limiting mutation rate and the evolutionary process underlying their extinction. Open questions include: for how many generations, and how large can populations grow at the limit mutation rate? What is the genetic structure of populations undergoing EEX? Can mutation drive the extinction of an exponentially growing population?

Here we propose a continuous-time multi-type critical birth–death process that allows theoretical exploration of the questions above, amongst other applications of biological interest. This models a population of cells that divide, die or mutate independent of each other, where the time between events is exponentially distributed. To mimic the behaviour of cells at the limiting mutation rate, we set the rates of division (α) and death (β) or mutation (ν) to be balanced (α=β+ν). In the simplest model for EEX, there is a single type of cells that divides (1)→(1)+(1) or acquires a lethal mutation (1)→∅ at the same rate. This is the case of haploid yeast and bacteria with mutation rates of one lethal mutation per cell division, which undergo EEX within a few generations (Morrison et al. 1993; Fijalkowska and Schaaper 1996; Herr et al. 2011). A more biologically interesting case is to allow multiple mutations to accumulate before a lethal mutation arrives. We model this by considering different cell types, where we call a cell type-i if it has accumulated i mutations. We assume that there is a maximal number n of mutations a cell can bear. Note that in this model mutations accumulate consecutively, therefore the lethal mutation is the nth mutation, not a particular mutation.

In the simplest multi-type critical process, each type-i cell divides independently (i)→(i)+(i) at rate one, or accumulates a new mutation (i)→(i+1) also at rate one, except type-n cells, which divide or die at rate one. We represent this birth-mutation process by the following scheme1

A more general version comes from introducing death and allowing different types of cells to divide and mutate at different rates. We refer to this as the n-type birth–death process, which can be illustrated as2

Note that for any given type, cells appear (via division) and disappear (via mutation or death) at the same rate, and thus all types remain critical. The overall process is critical, hence eventually all cells go extinct with probability one. This can be seen intuitively, since critical type-1 cells go extinct without a source, at which point the same can be said about type-2 cells, and so on. Figure 1 shows the evolution of the number of cells of each type in a four-type process of (1). Note that each type grows faster than the previous and takes longer to become extinct. Thus, the more mutations that can accumulate before a lethal mutation, the longer it takes for the population to undergo EEX. Interestingly, the population grows approximately exponentially during a transient, but eventually goes extinct with probability one. This behaviour cannot be captured by a supercritical process with an accumulation of deleterious mutations. In that case, if the initial types of cells have higher fitness, these will overtake any subsequent less fit types, resulting in a positive probability of survival of the population. Moreover, the n-type critical process does not require a reduction in fitness until the last lethal mutation, naturally capturing the phenomenon of EEX.Fig. 1 An example simulation of the process of (1) run until extinction with n = 4 types, where each coloured area represents the number of cells of a type, with different types piled on top of each other. Hence the envelope is the total population size. A single initial type-1 cell resulted in a short-lived type-1 population which seeded a type-2 population before its extinction. Note how type-3 and type-4 cells were initiated multiple times. One can observe the initial fast growth of the total population before extinction

Due to their applicability to model exponentially growing populations, multi-type, super-critical processes with consecutive mutations (decomposable) have been extensively studied (Athreya and Ney 2004; Kesten and Stigum 1967; Durrett 2015; Nicholson et al. 2022). In the critical case, the asymptotic behaviour has been studied by many authors (Sevast’yanov 1959; Chistyakov 1959; Mullikin 1963). Foster and Ney (1976) derived the asymptotic survival probability for discrete-time models. Similar results were obtained simultaneously by Ogura (1975), also including continuous time processes. Later on, Foster and Ney (1978) proposed limit theorems for the generating functions of population sizes, conditioned on the survival of the first type of cells, a condition less relevant for biological applications.

In this work, we derive large-time asymptotic solutions of multi-type critical birth–death processes exploiting the fact that, after an initial stage of exponential growth, the population is dominated by the last type. More precisely, we show that, for a large time, the system is non-empty with probability proportional to (time)-χn, where χn=21-n, establishing the relationship between the time-dependent survival probability and the maximal number of mutations that can accumulate. The exponents of the survival probability were derived by Foster and Ney (1976), Ogura (1975) through a different approach. Knowing the asymptotic survival probability allows us to find an appropriate scaling of the system, and derive the limiting distribution for the number of cells of a given type present at time t. We find asymptotic solutions for the distribution of cells of type k=1,⋯,n, as well as the total number distribution. Our methods extend the exact solution and limit results (Antal and Krapivsky 2011) for the two-type critical birth–death case, and complement the results of recent work on super-critical processes (Nicholson et al. 2022). Remarkably, we find that the distributions for the number of cells of a given type and the total number distributions have algebraic and stationary tails, described by the same exponents χn as the survival probabilities. This provides interesting biological insight into the behaviour of the modeled populations, showing that they can reach a stationary growth phase before extinction. We derive further estimates of interest for studying cancer and bacteria growth, including the distributions of time of arrival and extinction of cells with an arbitrary number of mutations.

The next sections progressively build up to the main results of this work, which can be found in Sect. 5. We introduce the simplest critical processes in Sects. 2 and  3. We present exact solutions to the two-type critical process in Sect. 4 and derive asymptotic solutions to the n-type case in Sect. 5. In Appendix E, we examine the more general processes including death and arbitrary birth and mutation rates, and show that the behavior is essentially the same as in the simplest critical birth–death process. Hence in the bulk of the paper, we limit ourselves to the simplest version. Example applications of the results for studying EEX in microbes and tumours, as well as connections to experimental work, can be found in Sect. 6.

Single type

In the simplest case, EEX occurs because cells acquire a lethal mutation at the same rate as cell division (Morrison et al. 1993; Fijalkowska and Schaaper 1996; Herr et al. 2011). This is represented by a single type critical process. Let us recall its basic properties (Athreya and Ney 2004). We denote by Zi(t) the number of type-i cells at time t. We start with a single type-1 cell and study Z1(t) via the generating functionZ1,1(x1,t)=E(x1Z1(t)|Z1(0)=1)=∑a≥0Pa(t)x1a,

where the first index of Z1,1 refers to the type of the initial cell while the second to the number of types considered, and Pa(t)=P(Z1(t)=a|Z1(0)=1) is the probability of having a cells at time t. This generating function satisfies the backward Kolmogorov equation ∂tZ1,1=(1-Z1,1)2 with Z1,1(x1,0)=x1, from which3 Z1,1(x1,t)=t(1-x1)+x1t(1-x1)+1.

By expanding the generating function around x1=0 we obtain the probability to have a cells at time t,4 Pa(t)=1a!∂x1aZ1,1(0,t)=1(1+t)2t1+ta-1

for a≥1, and the survival probability5 S1,1(t)=1-Z1,1(0,t)=1-P0(t)=11+t.

Hence the number of cells conditioned on survival has a geometric distributionZ1(t)|{Z1(t)>0}∼Geo11+t.

Notice that the average number of cells remains constant,EZ1(t)=∂x1Z1,1(1,t)=∑a≥1aPa(t)=1,

throughout the evolution. Although the probability of extinction for t→∞ is one, it takes a long time, since for T=inf{t:Z1(t)=0},ET=∫0∞P(T>t)dt=∫0∞S1,1(t)dt=∞.

It is well known (Athreya and Ney 2004) that conditioned on survival Z1(t)/t converges in distribution to an exponential random variable6 Z1(t)t|{Z1(t)>0}→Y1∼Expo(1)

for t→∞. This is immediate from the properties of a geometric distribution when taking the limit t,a→∞ limit with a/t=y constantPZ1(t)t>y|Z1(t)>0=tt+1ty→e-y=P(Y1>y).

One can also derive this from the generating function. First note that P(Z1(t)=a|Z1(t)>0)=Pa(t)/(1-P0(t)) for a≥1, hence the conditional generating function7 E(x1Z1(t)|Z1(t)>0)=∑a≥1Pa(t)1-P0(t)x1a=Z1,1(x1,t)-P0(t)1-P0(t).

The large a limit is related to x1≈1 so, in order to get a non-trivial limit, we write x1=1-p/t with p constant and notice that for a,t→∞, the right-hand side of (7) convergesZ1,1(1-p/t,t)-P0(t)1-P0(t)=1-p/t1+p→11+p.

The convergence of (7) requires that tPa(t)/(1-P0(t))∼t2Pa(t)→fY1(y) in order to get∑a≥1Pa(t)1-P0(t)x1a=∑a≥1tPa(t)1-P0(t)1-ptyt1t→∫0∞fY1(y)e-pydy=Ee-pY1

to converge to the Riemann integral, which is the Laplace transform of the density fY(y). The two limits are of course the same, hence8 limt→∞Z1,1(1-p/t,t)-P0(t)1-P0(t)=Ee-pY1=11+p.

By inverting (8), the Laplace transform of the density fY1(y), we indeed obtain that Y1∼Expo(1).

Infinite types

An interesting special case is obtained when we consider infinitely many types (n=∞) in the process described by scheme (1). Biologically, this represents a population at a limiting mutation rate (hence every sub-population goes extinct) but in which infinitely many mutations can accumulate (hence the total population grows exponentially).

Since n=∞, there is no extinction, but up to type m the process is identical to the m-type model. What becomes simple though in the infinite type model is the total number of cells: if we disregard the types of cells then Z(t)=∑i≥1Zi(t) is a Yule process with rate one (Athreya and Ney 2004). The generating function of a Yule process Z(x,t)=E(xZ(t)|Z(0)=1) can be obtained from the backward Kolmogorov equation ∂tZ=Z(Z-1) with Z(x,0)=x, which leads toZ(x,t)=xx+(1-x)et.

By expanding this around x=0, we obtain the mass function of the total number of cells in the infinite type process with a single initial type-1 cellΠs(t)=P(Z(t)=s|Zi(t)=δi1)=1s!∂xsZ(0,t)=e-t(1-e-t)s-1,

The mean number of cells grows exponentially as EZ(t)=∂xZ(1,t)=et. For finite multi-type critical branching processes, Πs(t) is not known in general. This is addressed in the next sections.

Two types

Consider now the simplest two-type critical branching process, which models a population in which a maximal of two mutations can accumulate. This could be interpreted as a population of diploid cells with a single locus for a lethal mutation, in which double allelic mutation is required to acquire the lethal phenotype. The two-type critical birth–death process is represented by9

where all steps occur at equal rates; we set these rates to unity for simplicity. The general case with death is considered Appendix A. While the single type birth–death process is easily soluble, the two-type birth–death process reduces to generally unsolvable Riccati equations, which were recently shown to be solvable for certain birth–death models with general rates (Antal and Krapivsky 2010, 2011). For the critical case, the situation slightly simplifies and the results are more explicit. Although the analytic solution for the two-type critical case was published in Antal and Krapivsky (2011), here we re-derive the solution in more detail. These exact results will be useful for checking the validity of the asymptotic results we derive in Sect. 5 for general n, which is the only available approach for n≥3.

Since the process is decomposable, the generating function for type-1 cells is the same as in the single type process given in (3). The generating function for both cell types, starting with a single initial type-i cell is10 Zi,2(x1,x2,t)=E(x1Z1(t)x2Z2(t)|Zj(t)=δij)

for i=1,2, and hence the initial conditions are 11a Z1,2(x1,x2,0)=x1

11b Z2,2(x1,x2,0)=x2.

A convenient setting for studying the two-type branching process is provided by the backward Kolmogorov equations 12a ∂tZ1,2=Z1,22+Z2,2-2Z1,2

12b ∂tZ2,2=Z2,22+1-2Z2,2.

Let σ be an operator that increases the indices of a function by one,13 σ(Fi,j(xi,xj))=Fi+1,j+1(xi+1,xj+1).

Note that in the general case (see Appendix E) σ increases the indices of the rates αi,χi,νi too. Since we are considering decomposable processes, one can easily see that Zi+1,j+1=σ(Zi,j), namely the j-type process starting with a single i cell is the same as the j+1 process starting with a single i+1 cell, up to a change in indices. In particular, Z2,2=σ(Z1,1), where Z1,1 is the generating function of the single type process, given in (3).

Survival probabilities

Before solving Eqs. (12a)–(12b), we consider the simpler task of finding the survival probability of the system, where ‘survival’ refers to the situation when there are alive cells of any type at time t,Si,2(t)=P(Z1(t)+Z2(t)>0|Zk(0)=δi,k).

The survival probabilities are related to the corresponding generating functions for type-2 cells, S1,2(t)=1-Z1,2(0,0,t) and S2,2(t)=1-Z2,2(0,0,t). Hence from Eqs. (12a)–(12b) we deduce that the survival probabilities S1,2(t) and S2,2(t) evolve according to rate equations 14a dS1,2dt=S2,2-S1,22

14b dS2,2dt=-S2,22,

with initial conditions S1,2(0)=S2,2(0)=1. Note that15 S2,2=σ(S1,1)=S1,1=11+t

is the survival of the single type process, given by (5). The governing Eq. (14a) therefore becomes an initial value problem16 dS1,2dt=11+t-S1,22,S1,2(0)=1.

This is a soluble Riccati equation. We first simplify the non-homogeneous term by using the new time variable τ=1+t, which leads todS1,2dτ=2τ-2τS1,22.

We transform this into a second-order linear differential equation by setting S1,2(τ)=12τB′(τ)B(τ) to getB′′(τ)-B′(τ)τ-4B(τ)=0.

Finally, we let s=2τ and write B(τ)=s2A(s) to obtain the standard differential equation for modified Bessel functions (NIST 2023)s2A′′(s)+sA′(s)-(1+s2)A(s)=0.

Hence the solution is a linear combination of these functionsA(s)=c1I1(s)+c2K1(s),

where Ip and Kp are the modified Bessel functions of order p of the first and second kind, respectively. Making the previous substitutions backward, we arrive at the solution of (16)17 S1,2=1τI0(2τ)-cK0(2τ)I1(2τ)+cK1(2τ),τ=1+t.

The initial condition S1,2(t=0)=1 fixes the amplitude18 c=I0(2)-I1(2)K0(2)+K1(2).

Generating functions

We now return to the generating functions Z1,2 and Z2,2. A full solution of Eqs. (12a)–(12b) subject to the initial conditions (11a)–(11b) is obtained in similar way as the above solution of Eqs. (12a)–(12b). First, notice that the solution of (12b) is just the generating function (3) of a single type but with the indexes increased by one19 Z2,2=σ(Z1,1)=t(1-x2)+x2t(1-x2)+1.

Plugging this into (12a) we recast it into a Riccati equation20 d(1-Z1,2)dt=1(1-x2)-1+t-(1-Z1,2)2.

Comparing (20) and (16) we see that 1-Z1,2 satisfies the same equation as S1,2, the only distinction is in the shift of the time variable in the first term on the right-hand sides. Thus we use the variable21 τ=t+(1-x2)-1.

and then the solution becomes22 Z1,2=1-1τI0(2τ)-cK0(2τ)I1(2τ)+cK1(2τ).

The initial condition Z1,2(x1,x2,0)=x1 fixes the amplitude23 c=I0(2τ0)-(1-x1)τ0I1(2τ0)K0(2τ0)+(1-x1)τ0K1(2τ0),τ0=(1-x2)-1/2.

For x1=x2=0 we recover the survival probabilities (17).

Having a full solution (22)–(23) for the generating function, we can extract the probability of finding a type-1 cells and b type-2 cells, Pa,b(t), by using Cauchy’s integral formula:24 Pa,b(t)=1(2πi)2∮dx1x1a+1∮dx2x2b+1Z1,2(x1,x2,t).

Unfortunately, explicit solutions for Pa,b(t) are not available for arbitrary a and b, although one can reduce the number of integrations in Eq. (24) to one, and use numerical methods to obtain values for arbitrary a and b, see Appendix B and C. Fortunately, in many applications, a less detailed description suffices. For instance, one may be interested in the marginal distribution of type-2 cells,Ps(2)(t)=P(Z2(t)=s)

or the total population size of the two-type processΠs(2)(t)=P(Z1(t)+Z2(t)=s).

These are encoded in the generating functions25 Z1,2(1,x,t)=E(xZ2(t)|Zj(0)=δ1j)=∑s≥0Ps(2)(t)xs

26 Z1,2(x,x,t)=E(xZ1(t)+Z2(t)|Zj(0)=δ1j)=∑s≥0Πs(2)(t)xs.

In the next section, we derive the asymptotic behaviour of both distributions.

n types

In this section, we study the n-type critical birth–death process. This models a population of cells that mutate at the limit mutation rate until they reach a maximal number of mutations n. This is represented by the scheme

This considers a decomposable critical process in which each type can only be produced by the preceding type, and all types divide or mutate at rate one, except the maximal type-n cells, which cannot mutate but die at rate one. The results are easily generalized to consider more general cases including death and multiple mutation rates; the derivations can be found in the Appendix E.

The naïve approach to the problem is to study EZi(t). The forward equations are ddtEZ1(t)=0 andddtEZi(t)=EZi-1(t)

for i=2,⋯n, with initial condition EZi(0)=δi1. The solution is27 EZi(t)=ti-1(i-1)!.

For the infinite type case (n→∞) we recover the Yule process of Sect. 3 for the total number of cells EZ(t)=∑i≥1EZi(t)=et. For a finite number of types, Z(t)=∑i=1nZi(t), and the mean total number of cells grows exponentially EZ(t)≈et for small t but algebraically as tn-1 when t→∞. Thus, even though the process is critical and dies out at a finite time with probability one, EZ(t) keeps growing forever, and therefore the naïve approach fails to capture one important aspect of the system. Hence instead of working with the mean number of cells, we derive asymptotic solutions for large time limits following a similar strategy introduced for the single type process.

The generating function for the n-type cell process starting with a single initial type-i cell is28 Zi,n(x1,x2,⋯xn,t)=E∏j=1nxjZj(t)|Zk(t)=δik

which can be obtained by solving the backward Kolmogorov equations29 ∂tZi,n=Zi,n2+Zi+1,n-2Zi,nfori=1,⋯,n-1Zi,n2+1-2Zi,nfori=n.

Thus, in order to obtain the interesting generating function Z1,n for the n-type process starting from a single type-1 cell, we need to first solve the n-1 equations for Zn,n,Zn-1,n,⋯,Z2,n. However, since the system is decomposable, we have that Zi+1,j+1=σ(Zi,j), where σ() is the index increase operator (13). Thus, if we know the generating function for type-k cells, solving the system for k+1 involves solving a single additional equation. For example, for the 3-type process, Z3,3=σ(Z2,2) and Z2,3=σ(Z1,2), which we solved in the 2-type process, and thus we only need to solve one more equation in order to obtain Z1,3. The same applies to the survival probabilities, which we consider in the following section.

Survival probabilities

As seen in Sect. 4, it is simpler to solve the system for the survival probabilities30 Si,n(t)=P(Z1(t)+⋯+Zn(t)>0|Zk(0)=δi,k)=1-Zi,n(0,0,⋯,0,t).

Substituting into (29), we deduce that31 dSi,ndt=Si+1,n-Si,n2fori=1,⋯,n-1,-Sn,n2fori=n,

with initial conditions Si,n(0)=1 for all i=1,⋯,n (see Fig. 2 for exact and numerical solutions). For any n, the solutions starting with a single type-n cell are given by the solutions of the single-type system up to shift in indices,Sn,n=σ(Sn-1,n-1)=⋯=σ(S2,2)=σ(S1,1)=1/(1+t).

To get Sn-1,n, we need to solve32 dSn-1,ndt=11+t-Sn-1,n2.

We have solved this case in Sect. 4.1 exactly, where we specified n=2, but the solution’s form is the same for all n. Even though analytic solutions are invaluable, the most interesting long-time asymptotic behavior can be extracted directly from the above equation using standard methods (Bender and Orszag 2013), circumventing the exact solutions. Assume that asymptotically Sn-1,n∼(1+t)-α, with α>0, where f(x)∼g(x) means f(x)/g(x)=1 for x→∞. Substituting into (32) we get-α(1+t)-α-1+⋯=(1+t)-1-(1+t)-2α+⋯

where the dots refer to smaller order terms. One needs to match the coefficients of the leading order terms on the two sides, which leads to a contradiction for both 2α>1 and 2α<1, thus 2α=1. That is, in the leading order, the right-hand side of (32) should vanish faster than 1/(1+t). This gives33 Sn-1,n∼(1+t)-1/2.

To get to higher orders we just add new terms to (33) of the form b(1+t)-β one by one with smaller values of β each time, substitute into (32) and match the coefficients. For the first 3 terms in the t→∞ asymptotic limit we get34 Sn-1,n∼(1+t)-1/2+14(1+t)-1+332(1+t)-3/2.

Here f(x)∼g(x)+h(x) means [f(x)-g(x)]/h(x)=1 for x→∞. Note that we performed the expansion in powers of 1+t for convenience, but one could simply expand the result in powers of t instead. This expansion can be also derived from the exact formula (17) using large argument asymptotic expansion of the Bessel functions (Appendix D).

Continuing to the next type, we have the differential equation for Sn-2,n in (31) where Sn-1,n is replaced by its above expansion. This procedure gives the leading order for all types by setting the time derivative to zero. In the leading order, we get that Sn-j,n∼(1+t)-χj+1 with χj=21-j. With a little more work, one obtains the next correctionSn-j,n∼(1+t)-χj+1+2-j-1(1+t)-12-χj+1

for j=1,⋯,n-1. The most interesting survival probability is S1,n, which gives the probability that the n-type system starting with a single type-1 cell is not empty at time t. Asymptotically, this is given by35 S1,n∼(1+t)-χn+χn2(1+t)-12-χn

for n≥2, and for n=1 the first term is exact. One can consider higher-order terms by successively adding terms and matching coefficients. Figure 3 shows that the second-order asymptotic accurately describes the behaviour.Fig. 2 Survival probabilities of the n-type critical process with initial condition Z1(0)=1. We plot the survival probability of the entire system, S1,n(t) (lines) and that of just type-n cells, Q1,n(t) (dashed-dots lines). Curves are for different types n=1,⋯,5. For type-1 and type-2 cells, the solutions are exact from (3) and (22), respectively. For other n>2, solutions are obtained numerically from (31) with different initial conditions

Fig. 3 Survival probability of the n-type process, for n=1,2,3,4,8. Second-order asymptotic solutions obtained from (35) (lines) are plotted together with results from simulations (dots)

Generating functions

We now attempt to find solutions for the generating functions for the number of cells. We rewrite the Kolmogorov equations (28) into36 d(1-Zi,n)dt=(1-Zi+1,n)-(1-Zi,n)2fori=1,⋯,n-1,-(1-Zi,n)2fori=n,

which for 1-Zj,n are identical to equations for the survival probability (31), but with initial conditions 1-Zj,n(x1,⋯,xn,0)=1-xj. This is not surprising since the survival probability is Si,n(t)=1-Zi,n(0,⋯0,t).

We again treat this system by first considering the generating function for the system starting with a single type-n cell, which corresponds to the single type case. Hence1-Zn,n=1t+(1-xn)-1,

which is the same as if we replace 1+t in Sn,n by t+(1-xn)-1. Starting with a single type-(n-1) cell we have37 d(1-Zn-1,n)dt=1t+(1-xn)-1-(1-Zn-1,n)2.

To get a nontrivial large-time behaviour we need (1-xn)-1∝t. Noting the similarity between (37) and (32), we assume a power behavior for 1-Zn-1,n which fixes the exponent and leads to1-Zn-1,n∼t+(1-xn)-1-1/2.

The same result can be also obtained directly from the explicit solution (22) for the generating function of type-2 cells Z1,2(1,x2,t) by taking the large argument asymptotic of the modified Bessel functions (Antal and Krapivsky 2011).

Continuing this procedure we get that, to the leading order,1-Zn-2,n∼t+(1-xn)-1-1/4,

and so on. In particular, the generating function starting with a single type-1 cell is given by38 1-Z1,n∼t+(1-xn)-1-χn.

Thus we see that the same exponents, χn=21-n, that describe the survival probabilities give the asymptotic behaviour of the generating functions. Note that in the leading order, the generating functions only depend on the last type of cells xn. This is not surprising if we note that the survival probability of type-k cells is the square root of the survival probability of type-k-1 cells. Thus the system becomes dominated by the last type. That is, the generating function for the last type is asymptotically the same as the generating function for the total number of cells,39 Z1,n(1,⋯,1,xn,t)∼Z1,n(xn,⋯,xn,t).

If we take (39) at xn=0 we see that the survival probability of just the type-n cells40 Qi,n(t)=P(Zn(t)>0|Zk(0)=δi,k)=1-Zi,n(1,⋯,1,0,t)

is asymptotically the same as the survival probability of S1,n(t) of the whole system, that is41 Q1,n(t)∼S1,n(t).

This can also be seen in Fig. 2, where we plotted exact (n=1,2) and numerical (n>2) solutions for S1,n(t) and Q1,n(t). It is clear from (29) that Qi,n(t) is governed by the same equations as the survival probability Si,n(t) of the entire system (31) but with initial conditions Qi,n(0)=0 for i=2,⋯,n-1 and Qn,n(0)=1. In Fig. 2 we see that subsequent types arise fast but disappear slowly, as the system gets dominated by the last type.

Cell number distributions

We now derive the asymptotic distribution of the number of type-n cellsPs(n)(t)=P(Zn(t)=s|Zk(0)=δk,1)

which is encoded in the generating function42 Z1,n(1,⋯,1,x,t)=E(xZn|Zk(0)=δk,1)=∑s≥0Ps(n)(t)xs.

To get non-trivial large-time behaviour, we need to condition the generating function on the survival of type-n cells. Following the procedure of Sect. 2, we express the conditional generating function as43 E(xnZn(t)|Zn(t)>0)=Z1,n(1,⋯,1,xn,t)-P0(n)(t)1-P0(n)(t).

As in the single type case, we obtain a nontrivial scaling by using the scaling variables xn=1-p/t and y=s/t when taking t→∞ with p and y constants. Substituting the leading order asymptotic expressions for the survival probability of type-n cells from (41) and (35)1-P0(n)(t)=Q1,n(t)∼S1,n(t)∼t-χn

and that of the generating function from (38)1-Z1,n(1,1-p/t,t)∼(t+t/p)-χn.

we get thatE(xnZn(t)|Zn(t)>0)→1-pp+1χn.

To obtain this convergence in a different way requires the following convergence to the density of a random variable YntPs(n)(t)Q1,n(t)∼t1+χnPyt(n)(t)→fYn(y)

so the generating function becomes a Riemann integral in the t→∞ limit:E(xnZn(t)|Zn(t)>0)=∑s≥1Ps(n)(t)Q1,n(t)xs=∑s≥1tPs(n)(t)Q1,n(t)(1-p/t)yt1t→∫0∞fYn(y)e-pydy=Ee-pYn=1-pp+1χn.

Hence we obtained a convergence in distributionZn(t)t|{Zn(t)>0}→Yn.

Inverting the above Laplace transform Ee-pYn we express the density of the limit variable via the confluent hypergeometric function (NIST 2023)44 fYn(y)=χnF(1+χn;2;-y).

Therefore, the scaling form of the type-n distribution is45 Ps(n)(t)≈χnt-1-χnF(1+χn;2;-s/t).

For n=1, since χ1=1 and F(2;2;-x)=e-x, we recover the exponential limit of the single type case given in (6).

Finally, for n≥2, by taking the large argument asymptotic of the confluent hypergeometric function (13.7.2 in NIST (2023)), we obtain the large y=s/t tail of the distribution.46 Ps(n)(t)≈χnΓ(1-χn)s-1-χnwhent≪s≪tn-11-χn.

The algebraic decay for n≥2 is in stark contrast with the exponential decay of the first type cells, n=1. Surprisingly, the tail is not only algebraic but also stationary. The algebraic tail describes the behaviour for a limited range of s values. The lower bound s≪t comes from the large y=s/t expansion. The upper bound is more subtle, and comes from noticing that the algebraic tail would imply an infinite mean, in contradiction with the finite mean we obtained in (27). To reconcile this, we set an upper bound, s∗, and estimate its order of magnitude,EZn∝tn-1∝∫1s∗ss-1-χnds∝s∗1-χn

which gives s∗∝tn-11-χn as we announced.

The validity of the scaling limit (45) is illustrated in Fig. 4 for n=2 via comparison to numerical solutions. One can see how the range of the algebraic tail expands with time. In Fig. 5, the scaled number distribution χn-1t1+χnPs(n)(t)=F(1+χn;2;-s/t) as given by (45) is compared to simulations as a function of s/t, where we chose t=20. For n>1 the asymptotic solution matches simulations, but eventually overestimates the probability for large s. For larger t, the asymptotic solution is better at the large s limit but under-estimates for small s. Within the derived range, the cell number distribution is well described by stationary tail (size)-1-χn.

We have derived the asymptotic behaviour of the last type of cells in the n-type process. Since the system is decomposable, in order to get the distribution of the previous types k=1,2,⋯,n-1, we simply need to stop the process at the kth type. In biological applications, we may also be interested in the total number of cells Z=Z1+⋯+Zn, with distributionΠs(n)(t)=P(Z(t)=s|Zk(0)=δ1,k).

This is encoded in the generating function47 Z1,n(x,⋯,x,x,t)=∑s≥0Πs(n)(t)xs,

hence we see by (39) that it is asymptotically the same as the distribution of the last type of cells,Πs(n)(t)∼Ps(n)(t).

Thus, the above analysis gives us access to both the asymptotic behaviour of the total number of cells, as well as the behaviour of each individual cell type.Fig. 4 Scaling of the probability Ps(2)(t) of finding s type-2 cells at time t; different lines indicate different times t. Probabilities are calculated numerically from the exact generating function (22) via the Inverse Fast Fourier Transform algorithm, see Appendix C. In the double limit t,s→∞, with s/t constant, the scaled distributions converge to the scaling limit given by (45), depicted by dashed line

Fig. 5 The scaled number distributions χn-1t1+χnPs(n)(t) from simulations (dots) are compared to the theoretical prediction F(1+χn;2;-s/t), cf. Eq. (45), as a function of s/t, where we fix t=20. Curves are for n=1,2,3,4. Note the power law decay for n≥2 as opposed to the exponential decay for the single type case (n=1)

Arrival and exit times

Let us study when new cell types appear and disappear from the system. It is easier to work with the infinite type version of the model for this question, otherwise, we need to assume that the type in question is not greater than n. For convenience, we study when type-n cells appear or disappear, but the results are of course valid for any type, not only the last.

For the pure birth-mutation process, the exit time of types is straightforward. LetEn=sup{t≥0:Zn(t)>0}

denote the time of extinction of type-n cells. For the birth-mutation process, notice that{En>t}={Z1(t)+⋯+Zn(t)>0}.

Hence the distribution of En when starting from a single type-1 cell is given by the total survival of the n-type process, which we derived asymptotically in (35).

Let us now turn to the arrival timeTn=inf{t≥0:Zn(t)>0}

when the first type-n cell appears. We are interested in its distribution starting with a single type-i cellhi,n(t)=P(Tn>t|Zk(0)=δik).

To derive an equation for this quantity consider a modified system where type-n cells neither divide nor die, just stay alive forever. Hence their generating function, when starting with a single type-n cell, stays constant Zn,n=xn, but all other equations for i<n remain the same in (29). The existence of a type-n cell in the modified system then indicates that they were produced also in the original system. Hence hi,n(t)=Zi,n(1,⋯,1,0,t), and thus48 dhi,ndt=hi,n2-2hi,n+hi+1,n

with initial condition hi,n(0)=1 for i=1⋯,n-1 and hn,n(t)≡0. In terms of gi,n(t)=1-hi,n(t)=P(Tn≤t|Zk(t)=δik), Eq. (48) takes a simpler form49 dgi,ndt=-gi,n2+gi+1,n

with initial condition gi,n(0)=0 for i<n and gn,n(t)≡1 for all t.

For the arrival of the first type-2 cell, the above equation becomes50 dg1,2dt=-g1,22+1

with solution51 g1,2(t)=tanht.

Hence the first type-2 cell arrives on average atET2=∫0∞1-g1,2(t)dt=log2.

Note also that52 h1,2(t)=1-g1,2(t)=P(T2>t|Zk(0)=δ1k)=21+e2t

have the same form as for the analogous supercritical process (Nicholson et al. 2022).

Now we can use this solution for the arrival of type-3 cells, by noting that gi+1,n=σ(gi,n-1) and we get53 dg1,3dt=-g1,32+tanht.

The solution of (53) can be expressed as a lengthy combination of hypergeometric functions. For subsequent types, no exact solutions are available.

Let us study instead the more general birth–death process described by scheme (2). The above method presented for the birth-mutation process stays valid, and for all αi=1 and a constant mutation rate ν for all types, Eq. (49) becomes54 dgi,ndt=-gi,n2+νgi+1,n

with initial conditions gi,n(0)=0 and gn,n≡1. The simplest non-trivial property of this process is the probability that the first type-n mutant eventually arrives starting from a single type-1 cell,g1,n(∞):=limt→∞g1,n(t).

For the simple birth-mutation process (ν=1) this quantity is trivial: g1,n(∞)=1 for all n, namely all types arrive eventually with probability 1. This is not the case in the more general birth–death case where, by setting the left-hand side of (54) to zero, we obtain thatg1,n(∞)=ν1-χn.

If we denote by M the maximal type that ever appears, then what we found is that P(M≥n)=g1,n(∞)=ν1-χn. It gets easier for larger types to appear, in the sense that P(M≥n+1|M≥n)=ν-n→1 as n→∞. This also implies that the mean number of types that ever appear is infinite, EM=∞.

Since g2,2≡1, we can get an explicit solution for the arrival of the first type-2 cellg1,2(t)=νtanhνt.

which generalises (51).

For the arrival of the first type-3 cell, using that g2,3=σ(g1,2)=g1,2, we have that55 dg1,3dt=-g1,32+ννtanhνt.

The initial condition is g1,3(0)=0. Let us normalize g1,3 by its limiting probability g1,3(∞)=ν3/4. That is, we introduce g~1,3:=g1,3g1,3(∞) to get1g1,3(∞)dg~1,3dt=-g~1,32+tanhνt.

Next, we re-scale time, t→t~=tg1,3(∞)=tν3/4, and definek1,3(t~)=g~1,3t~ν-3/4=ν-3/4g1,3(t)

to arrive atdk1,3dt~=-k1,32+tanhν-1/4t~.

In the limit where ν→0 we getdk1,3dt~≈-k1,32+1,

which is the same differential equation we have for g1,2 for ν=1 in (50) with the same initial condition k1,3(0)=0. Hence k(t~)=tanht~ and theng1,3(t)≈ν3/4tanhν3/4t.

The same procedure can be applied to obtain the arrival distribution of all types,56 g1,nt≈g1,n(∞)tanh(g1,n(∞)t)=ν1-χntanhν1-χnt

which can be then verified by induction. In Fig. 6 we observe that, for a fixed type, the asymptotic agrees with the behaviour in the ν→0 limit. However, for a fixed ν, the approximation becomes worse as we consider more types.Fig. 6 On the left, the scaled probability of arrival of the first type-3 cell, ν-3/4P(T3≤t) in terms of ν3/4t, obtained numerically from (54) for different ν is compared to the scaling limit (56). On the right, the normalized arrival probability νχn-1P(Tn≤t) of different types obtained numerically from (54) are compared to the corresponding scaling limit (56) (dashed line) for ν=10-5 and n=2,⋯,5

Examples and applications

We now discuss example applications of the results for studying the evolution of microbes and tumour cells at the limiting mutation rate and relate them to experimental work.

Colony size and number of generations until EEX in microbes

In bacteria, haploid and diploid yeast, numerous experimental projects have investigated EEX by crossing cell lines with different mutator alleles Morrison et al. (1993), Fijalkowska and Schaaper (1996), Herr et al. (2011), Herr et al. (2014), Soriano et al. (2021). The observable quantities of such experiments are the CanR mutation rate of the cell lines (number of mutations in CanR gene per division), and the number of viable colonies after a given number of generations. The main goal is to identify the maximum mutation rate. We have seen that defining the limiting mutation rate is conceptually straightforward: if the sum of death and mutation rate is at least equal to the division rate, and there is a finite number of mutations that can accumulate, then EEX will occur with probability one. However, for how long and how large the colonies can grow depends on a number of factors, including the number of mutations that can accumulate before the lethal mutation. These are questions of interest for experiment design, which we can answer using the model proposed.

If we set the birth rate to 1 in the no-death model, we can interpret time in units of cell divisions or generations. We can immediately compute the probability that a colony initiated by a single cell will survive after t generations immediately from Eq. (35), which gives the survival probability after t generations depending on the maximal number of mutations that can accumulate n,S1,n∼(1+t)-χn+χn2(1+t)-12-χn

where χn=21-n. In Fig. 7a, we see that for n=1, the survival probability is less than 10% after 10 generations, whilst up to 100 generations are needed for the same probability if n=2. For n=5, a reduction in survival probability of 50% will take 105 generations. In yeast, most experiments determine the growth of colonies after 20 generations. We can compute the expected size of the colony after t generations using Eq. (47). In Fig. 7b, we see that, if only one mutation can accumulate (n=1), we won’t observe any growth after 20 generations, whilst if two or three mutations can accumulate (n=2,3), colonies will form, although smaller than expected. This is qualitatively relatable to the comparison between haploid and diploid yeast investigated by Herr et al. (2014), who found that at the limiting mutation rate of one mutation in an essential gene per cell division, haploid yeast colonies are not viable, whilst diploid yeast form viable, but smaller colonies. For n>4, the expected growth after 20 generations plateaus to the expected size in a normally exponentially growing population, suggesting that one would need to run longer experiments to observe EEX at the macroscopic level in cell lines that are more robust to mutational burden, in accordance with 7A. The survival probabilities and expected population size for the general case including cell death are available in Appendix E). Introducing death in the model effectively re-scales time, such that sub-populations with higher death rates go extinct faster.Fig. 7 (a) Number of generations until the survival probability is less than 10% and 50%, as a function of the maximal number of mutations that can accumulate, calculated from (35). (b) The expected size of colonies started with one single cell after 20 generations, as a function of the maximal number of mutations that can accumulate, calculated from (47). The dashed grey line shows the expected size of a population growing exponentially with no maximal number of mutations

A remarkable conclusion from the model is that, even if cells are at the limiting mutation rate and hence EEX will occur with probability one, this might take many generations. Thus, in the usual experimental setting, EEX is not detectable by following colony size. However, it would be detectable if sequencing or genetic information of the colonies at different time points was available. In the next Section, we consider the evolution of genetic diversity in more depth.

Genetic diversity during EEX

Multi-type branching processes are widely applied to studying the genetic structure of growing populations. In the multi-type critical process, we observe an interesting behaviour, where extreme mutation rates result in an initial increase of genetic diversity, but after a transient phase, the population becomes dominated by cells carrying the maximal number of mutations, and thus loses genetic diversity. In Fig. 8 we plot the evolution of the number of cells of different simulation runs of the 3, 4, and 5-type process, next to the Shannon diversity index, which is defined as57 H′(t)=-∑i=1npi(t)logpi(t)

where pi(t) is the proportion of cells of type-i at time t. Indeed, we can see that H′ increases as the first mutations accumulate, and rapidly declines when the last type of cells arrive.Fig. 8 Example simulation runs of the 3-type (a), 4-type (b, c) and 5-type (d) process, where each colored area represents the number of cells of a type, with different types piled on top of each other, and we ran the process until extinction or 104 steps. Next to each run, we plot the corresponding Shannon diversity index H′(t) in time 57. In all examples, this is maximal when the last type arrives, and then quickly decreases as the population becomes dominated by the last type

Even though an analytical expression for H′(t) from the model is not available, we can use (35) to obtain the expected total number of types or subclones present at time t, with its large time asymptotic expression58 Kn(t):=E[#of types at timet]=∑inQ1,i(t)∼∑inS1,i(t)∼t-χn,

where Q1,i(t) denotes the survival probability of the ith cell-type, and S1,i(t) denotes the survival probability of the entire population (defined in (40) and (30), respectively). In Fig. 9 we see that, similarly to the Shannon diversity index, Kn(t) initially increases and then quickly decreases. As observed in the simulations, after 10–20 generations, the expected number of types present already declines, in sharp contrast with the case in which the colony grows exponentially at a high mutation rate but cells do not acquire lethal mutations. Therefore, even though the macroscopic behaviour is indistinguishable in that timescale, if genetic information on the colonies in time is available experimentally (e.g. the number of subclones after 20 generations), one can identify populations that will eventually reach extinction due to high mutation. Moreover, one can use the expressions derived from the model to infer the maximal number of mutations n, which can then be used to quantitatively predict the evolution of the colonies and the time until EEX.Fig. 9 Mean number of types present at time t for populations with different maximal number of mutations (n), calculated from (58) (solid lines) with their large time asymptotic expressions (dashed lines). The dotted grey line shows the expected number of types in a population growing exponentially with no maximal number of mutations (n=∞)

Discussion

Multi-type branching processes provide a natural tool to model biological processes driven by cell division, death, and mutation. Due to their potential to describe evolutionary dynamics, extensive work has been dedicated to deriving solutions of multi-type processes, especially the super-critical and sub-critical cases (Durrett 2015; Nicholson et al. 2022). In this work, we have focused on finite-type critical processes, which mimic populations in which mutations accumulate until a maximal number of alterations is reached, resulting in extinction.

Driven by biological motivation, we have focused on deriving solutions for the survival of the population, the number of cells, and the arrival and extinction time of cells with different mutations. We found that the survival probability of the overall system, which is asymptotically equivalent to the survival probability of the last type of cells, decays as t-χn for the nth type. The χn=21-n exponents for the survival probability had been derived by previous work by Foster and Ney (1976) and Ogura (1975), although following different approaches. With our approach we have also derived higher-order terms, facilitating the estimation of the accuracy of the first-order term.

For two cell types, the generating function of the population sizes was expressed explicitly in terms of modified Bessel functions. This provides a numerically efficient way to extract the number distributions, as detailed in Appendix C. By conditioning on survival of the final type, we extract the distribution of the number of type-k cells in the large time limit. The survival probabilities show us that the system becomes quickly dominated by the last type of cells, and indeed the distribution of the last type coincides with that of the total number of cells. This distribution only depends on the ratio s/t of the size s of the population and time t. Interestingly, in the large s/t limit, this has algebraic and stationary tail (size)-1-χn. That is, for a fixed (large) number of cells, there is a time regime during which the probability of finding s cells remains constant in time. Since this fat tail – for each but the first cell type – would imply infinite population size, we have derived an upper cutoff for the power law tail which ensures the finiteness of the mean population sizes. These power-law tails appearing for large times for n≥2 cell types are in sharp contrast to the purely exponential behaviour of the mass function of the first cell type. Although the solution is only valid for large time and population size, this corresponds to the range of interest in biological applications.

We provide exact or asymptotic formulas for our model of EEX in time, including the evolution of population size, and the time of arrival and extinction of sub-populations. In Sect. 6, we discussed example applications of these formulas for studying EEX in microbes, in relation to experimental work. In particular, we showed that populations at the limit mutation rate that can accumulate a high number of mutations take longer to exhibit macroscopic changes in colony size than the number of generations that experiments usually track. Thus, even though such cell lines would eventually undergo error-induced extinction, this extreme fate could not be predicted by analyzing colony growth. In order to identify colonies at the limit mutation rate, experiments should complement macroscopic measures with measures of genetic diversity such as the number of subclones.

Apart from providing tools to quantitatively analyze and guide experimental design, the model provides theoretical insight into error-induced extinction. We propose a general mechanism for EEX: populations of cells that divide and mutate or die at the same rates and have a maximal number of mutations tolerable. In the context of cancer, this maximal number might represent the amount of DNA damage that can accumulate before being detected by the immune system (Schumacher and Schreiber 2015). We find that the modeled populations become dominated by cells carrying the maximal number of mutations and thus lose genetic diversity. This is related to the idea that genetic instability results in extinction because populations cannot overcome selective barriers (Tejero et al. 2016; Andor et al. 2017; Tilk et al. 2022). Another interesting behaviour is that populations at the error threshold reach a stationary phase before going extinct. Stationary growth of cancer or bacteria populations is normally associated with having reached a carrying capacity (due to limited nutrients, etc.) (Gerlee 2013; Tjørve and Tjørve 2017). Our model shows that this macroscopic behaviour might be caused by mutational burden, in which case the fate of the population is drastically different, highlighting the importance of considering genetic structure when modeling population growth.

A related phenomenon to EEX is the so-called error catastrophe, which describes the inability of a genetic element to be maintained in a population as the fidelity of its replication machinery decreases beyond a certain threshold value, such that it cannot produce enough copies of itself (Summers and Litwin 2006). This has been invoked as a theoretical basis for the treatment of viral infection with drugs that would push the error rate for copying the viral genome beyond this threshold (Vignuzzi et al. 2005). In our model, every type of cells eventually disappears due to high mutation, and thus undergoes an error catastrophe. Error catastrophes might result in the most fit cells to be replaced by lower fitness population, but, unlike EEX, not necessarily in extinction. This more general case can be studied using the proposed model, but either allowing for infinite number of mutations, or adjusting the birth and death rate of the last type of cells. As shown in previous work, the applications of multi-type critical processes extend beyond error catastrophe and EEX, e.g., to modeling infectious disease spread (Antal and Krapivsky 2012).

Throughout this paper, we have focused on the simplest n-type critical process, with zero death rate for all but the last cell type. The more general case including cell death, and arbitrary mutation and division rates can be found in Appendix E. This could be relevant for modelling systems in which mutations result in cells with different division, mutation, and death rates, as long as all cell types have critical growth. Including intermediate types with sub-critical and super-critical growth remains a challenge for future work. Future work should consider non-consecutive mutations, allowing for specific mutational paths to be modeled. This is particularly relevant to the study of the fitness of mutator alleles in cancer evolution, which is typically driven by the accumulation of specific mutations.

Appendix A. Two types with general rates

The solution for the two-type case has appeared in Antal and Krapivsky (2011), but we recall the results here for ease of reference with α1 explicitly included in the formulas59 Z1,2(x1,x2,t)=1-ν1α2τI0(2τ)-cK0(2τ)I1(2τ)+cK1(2τ)

where we used the variableτ=ν1α2α1α2(1-x2)+tα1

and the amplitude60 c=I0(2τ0)-α2ν1-1τ0(1-x1)I1(2τ0)K0(2τ0)+α2ν1-1τ0(1-x1)K1(2τ0)

where τ0≡τ(t=0).

Appendix B. Integral formulas Pa,b(t)

The variable x1 appears in (22) only through c, while c is a ratio of linear (in variable x1) functions. Therefore Z1,2 itself can be transformed into the ratio of linear in x2 functions. We can expand in x1 and effectively we do not need to take the integral in the (complex) x1 plane. Thus for a=0 we have61 P0,b(t)=12πi∮dx2x2b+1Z1(0,x2,t)

while for a>062 Pa,b(t)=12πi∮dx2x2b+1τ0aτ(f′)a-1(fg′-f′g)(f+τ0f′)a+1.

Here we used shorthand notationf=I1(2τ)K0(2τ0)-K1(2τ)I0(2τ0)f′=K1(2τ)I1(2τ0)-I1(2τ)K1(2τ0)g=I0(2τ)K0(2τ0)-K0(2τ)I0(2τ0)g′=K0(2τ)I1(2τ0)-I0(2τ)K1(2τ0).

Appendix C. Numerical solutions for Pa,b(t) and Ps(2)(t)

Explicit solutions for Pa,b(t) are not available. Numerically, however, one can access the probabilities. Using Mathematica’s SeriesCoefficient command we getPa,b(1)≈0.1691020.1903910.07890050.0341278⋯0.156590.05295940.02247250.00989238⋯0.06925590.02986490.01383360.00643762⋯0.03063020.01605770.008105750.00399019⋯⋮⋮⋮⋮

for a,b=0,1,2,3, which values can be used to check simulations. Simulating the process (for a few seconds on an average laptop with a code in C) over 108 runs we getPa,b(1)≈0.1691070.1904100.0788730.034104⋯0.1566390.0529870.0224640.009880⋯0.0692770.0298860.0138260.006427⋯0.0306050.0160390.0081020.003988⋯⋮⋮⋮⋮

However, this approach only works for small t. In order to obtain numerical solutions for large t, one can invert the generating function via Fourier transform, and apply the Inverse Fast Fourier Transform (IFFT) algorithm as an efficient method to calculate the probability density from the generating function (Abate and Whitt 1992). As an example, we illustrate how to obtain Ps(2)(t), that is, the probability of there being s type-2 cells. In order to obtain the joint distribution, one needs to invert Eq. (62). Recall the generating function for type-2 cellsZ1,2(1,x,t)=∑s≥0Ps(2)(t)xs.

The correspondence to the Fourier transform is easier to see if we consider z=1/x,Z1,2(1,1/z,t)=∑s≥0Ps(2)(t)z-s

this has inversion63 Ps(2)(t)=12πi∮Cdzzs-1Z1(1,1/z,t)

where the contour C goes counterclockwise around the origin in the complex y plane, and must enclose all poles of Z1(1,1/z,t). We choose C to be a circle of radius R enclosing all poles, and set z=Reiθ,64 Ps(2)(t)=Rs2π∫02πdθZ1(1,Re-iθ,t)eiθs.

We notice that we recover a Fourier series expansion. The above integral can be approximated as a sum by splitting the circle into N equidistant discrete points65 ps(2)(t)=RsN∑k=0N-1Z1(1,Re-2ikπ/N,t)e2iskπ/N

which is the Inverse Discrete Fourier Transform scaled by Rs. One can use the IFFT algorithm to obtain the approximate probability density function. The approximation error of ps(2)(t)≈Ps(2)(t) depends on the discretization (Abate and Whitt 1992; Cavers 1978). We use the following algorithm to extract the probabilities numerically: Calculate Z1(1,Re-iθ,t)=Z1(1,1/z,t),z=Reiθ.

Discretize the circle into N equidistant points.

Evaluate the generating function at each point. That is, calculate Z1[k]:=Z1(1,Re-2πk/N,t) for k=0,⋯,N-1.

Calculate the IDFT of Z1[k] using Fast Inverse Fourier Transform (e.g. fft.ifft in Numpy or InverseFourierTransform in Mathematica). This outputs N coefficients q0,⋯,qN-1.

Re-scale the sth coefficient by Rs to recover the probability of s cells, Ps(2)(t)≈ps(2)(t)=Rsqs.

Note that the generating function of the two-type system Z1,2(x1,x2,t) has a singularity at x2=1. Thus, we must take R>1. From Eq. (64) we see that the radius R does not affect the inversion. However, numerical errors due to finite precision number representation result in numerical errors with R. In Fig. 4 we take R=1.0001 as we find that if R is increased it becomes difficult to resolve the separate contributions from the poles and zeros of Z1(1,1/z,t). However, in other settings, numerical problems might arise if R is decreased to too close to the furthest pole. In general, it is common practice to take R to be 10% larger than the largest pole (Cavers 1978).

Appendix D. Useful formulas and asymptotic results

We use the large z expansions for the first and second type modified Bessel functions (Sect. 3.13 of Bender and Orszag (2013), Sect. 10.7 of NIST (2023)),66 In(z)∼ez2πz1-4n2-18z+(4n2-1)(4n2-9)2!(8z)2-⋯,

67 Kn(z)∼π2ze-z1-4n2-18z+(4n2-1)(4n2-9)2!(8z)2-⋯.

We often use I0(z)/I1(z)∼1, K1(z)∼0 and K0(z)∼0. On some occasions, we are interested in higher-order terms, which can be obtained by considering more terms in the asymptotic expansions. In particularI0(z)I1(z)=1+1z+12z2+2716z3+141128z4+Oz-5.

Appendix E. n types with general rates

Here we generalize our results to the n-type process to include death and mutation at arbitrary rates. This is described by the scheme

Let Zi,n(x1,x2,⋯xn,t) denote the generating function starting with a single initial type-i cell. The backward Kolmogorov equations read∂tZi,n=αiZi,n2+αi-νi+νiZi+1,n-2αiZi,nfori=1,⋯,n-1αn(Zn,n2+1-2Zn,n)fori=n.

E.1. Survival probabilities

We solve the system for the survival probabilities Si,n(t)=1-Zi,n(0,0,⋯,0,t), given bydSi,ndt=νiSi+1,n-αiSi,n2fori=1,⋯,n-1,-αnSn,n2fori=n,

with initial conditions Si,n(0)=1 for all i=1,⋯,n. Following the procedure outlined in Sect. 5.1 for the no-death case, we derive the survival probability of the n-type processS1,n(t)∼∏i=1n-1νiαiχi+1(αnt+1)-χn

where χn=21-n. Note that this system is equivalent to the multi-type critical process without death, up to a re-scaling of time. The asymptotic solutions are exact to the same order, but the convergence is slower.

E.2. Generating functions and total number distribution

We now turn to the generating functions. Again, we already know thatZn,n=1-1αnt+(1-xn)-1

and we know that for t→∞, in the leading order Zi,n only depends on the last type xn, and xn→1 so we may assume that in the leading order1-Zi,n∼(Bi(t+(1-xn)-1))-Ai

and match coefficients as usual to recover the coefficients Bi and Ai for all types. We arrive at the following expression for the leading order asymptotic generating function starting with a single type-1 cell,1-Z1,n(x1,⋯,xn,t)∼∏i=1n-1νiαiχi+1(αnt+(1-xn)-1)-χn.

Now we follow the procedure of the no-death case outlined in Sect. 5.2, use the scaling variables xn=1-p/αnt and y=s/αnt and take the limit t→∞ with p and y constants, to findE(xnZn(t)|Zn(t)>0)→∫0∞fYn(y)e-pydy=1-pp+1χn.

where we established the convergence in distributiontPs(n)(t)Q1,n(t)→fYn(y),Zn(t)t|{Zn(t)>0}→Yn

Hence we find that the limit random variable Yn does not depend on the parameters of the system, and has the density (44). Therefore, the scaling form of the n-type distribution is68 Ps(n)(t)∼χnQ1,n(t)tαnF(1+χn,2,-s/tαn)∼∏i=1n-1αiνiχi+1αntχn-1χnF(1+χn,2,-s/tαn).

Finally, by taking the large argument asymptotic of the confluent hypergeometric function (Appendix D), we obtain the stationary distribution69 Ps(n)(t)≈χnαnχn-1∏i=1n-1αiνiχi+1Γ(1-χn)s-1-χn.

As seen before, the algebraic stationary tail describes the behaviour up to an upper-cutoff s⋆ in the number of cells, which we recover by matching the order of magnitude of the mean EZn(t),70 t≪s≪tn-11-χn.

E.3. Arrival times

We now consider the arrival times for different birth and mutation rates for each type, αi, νi (as described by the scheme (2)). The equations for gi,j satisfy71 dgi,jdt=-αigi,j2+νigi+1,j

with initial conditions gi,j(0)=0 and gj,j≡1. This leads to the asymptotic probability of arrival of the jth type starting with a single cell72 g1,j(∞):=limt→∞g1,j(t)=∏i=1j-1νiαi-i.

As in the pure birth-mutation case (56), we obtain a small mutation rate limit (all νi→0, starting with νj-1) by normalizing and re-scaling time by (72)g1,jt≈g1,j(∞)tanhg1,j(∞)t.

Funding

Funding was provided by Centre for Doctoral Training in Mathematical Modelling, Analysis and Computation, University of Edinburgh (GB) (Grant No. EP/S023291/1).

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
==== Refs
References

Abate J Whitt W The fourier-series method for inverting transforms of probability distributions Queueing Syst 1992 10 5 87 10.1007/BF01158520
Abate J, Whitt W (1992) The fourier-series method for inverting transforms of probability distributions. Queueing Syst 10:5–8710.1007/BF01158520
Albertson TM Ogawa M Bugni JM Hays LE Chen Y Wang Y Treuting PM Heddle JA Goldsby RE Preston BD Dna polymerase ε and δ proofreading suppress discrete mutator and cancer phenotypes in mice Proc Natl Acad Sci 2009 106 40 17101 17104 10.1073/pnas.0907147106 19805137
Albertson TM, Ogawa M, Bugni JM, Hays LE, Chen Y, Wang Y, Treuting PM, Heddle JA, Goldsby RE, Preston BD (2009) Dna polymerase and proofreading suppress discrete mutator and cancer phenotypes in mice. Proc Natl Acad Sci 106(40):17101–1710419805137 10.1073/pnas.0907147106
Andor N Maley CC Ji HP Genomic instability in cancer: teetering on the limit of tolerance Cancer Res 2017 77 9 2179 2185 10.1158/0008-5472.CAN-16-1553 28432052
Andor N, Maley CC, Ji HP (2017) Genomic instability in cancer: teetering on the limit of tolerance. Cancer Res 77(9):2179–218528432052 10.1158/0008-5472.CAN-16-1553
Antal T Krapivsky PL Exact solution of a two-type branching process: clone size distribution in cell division kinetics J Stat Mech Theory Exp 2010 07 P07028
Antal T, Krapivsky PL (2010) Exact solution of a two-type branching process: clone size distribution in cell division kinetics. J Stat Mech Theory Exp 07:P07028
Antal T Krapivsky PL Exact solution of a two-type branching process: models of tumor progression J Stat Mech Theory Exp 2011 2011 08 P08018 10.1088/1742-5468/2011/08/P08018
Antal T, Krapivsky PL (2011) Exact solution of a two-type branching process: models of tumor progression. J Stat Mech Theory Exp 2011(08):P0801810.1088/1742-5468/2011/08/P08018
Antal T Krapivsky PL Outbreak size distributions in epidemics with multiple stages J Stat Mech Theory Exp 2012 07 P07018
Antal T, Krapivsky PL (2012) Outbreak size distributions in epidemics with multiple stages. J Stat Mech Theory Exp 07:P07018
Athreya KB Ney PE Branching processes 2004 New York Dover Publications
Athreya KB, Ney PE (2004) Branching processes. Dover Publications, New York
Bender CM Orszag SA Advanced mathematical methods for scientists and engineers I: asymptotic methods and perturbation theory 2013 Berlin Springer Science & Business Media
Bender CM, Orszag SA (2013) Advanced mathematical methods for scientists and engineers I: asymptotic methods and perturbation theory. Springer Science & Business Media, Berlin
Cavers JK On the fast fourier transform inversion of probability generating functions IMA J Appl Math 1978 22 3 275 282 10.1093/imamat/22.3.275
Cavers JK (1978) On the fast fourier transform inversion of probability generating functions. IMA J Appl Math 22(3):275–28210.1093/imamat/22.3.275
Chistyakov VP Generalization of a theorem for branching processes Theory Probab Appl 1959 4 1 103 106 10.1137/1104008
Chistyakov VP (1959) Generalization of a theorem for branching processes. Theory Probab Appl 4(1):103–10610.1137/1104008
Durrett R Branching process models of cancer 2015 Berlin Springer
Durrett R (2015) Branching process models of cancer. Stochastics in biological systems. Springer, Berlin
Fijalkowska IJ Schaaper RM Mutants in the exo i motif of escherichia coli dnaq: defective proofreading and inviability due to error catastrophe Proc Natl Acad Sci 1996 93 7 2856 2861 10.1073/pnas.93.7.2856 8610131
Fijalkowska IJ, Schaaper RM (1996) Mutants in the exo i motif of escherichia coli dnaq: defective proofreading and inviability due to error catastrophe. Proc Natl Acad Sci 93(7):2856–28618610131 10.1073/pnas.93.7.2856
Foster J Ney P Decomposable critical multi-type branching processes Sankhyā Indian J Stat Ser A 1976 38 1 28 37
Foster J, Ney P (1976) Decomposable critical multi-type branching processes Sankhyā. Indian J Stat Ser A 38(1):28–37
Foster J Ney P Limit laws for decomposable critical branching processes Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete 1978 46 1 13 43 10.1007/BF00535685
Foster J, Ney P (1978) Limit laws for decomposable critical branching processes. Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete 46(1):13–4310.1007/BF00535685
Fox EJ, Loeb LA (2010) Lethal mutagenesis: targeting the mutator phenotype in cancer. In: Seminars in cancer biology, Elsevier, vol 20, pp 353–359
Gerlee P The model muddle: in search of tumor growth laws Cancer Res 2013 73 8 2407 2411 10.1158/0008-5472.CAN-12-4355 23393201
Gerlee P (2013) The model muddle: in search of tumor growth laws. Cancer Res 73(8):2407–241123393201 10.1158/0008-5472.CAN-12-4355
Herr AJ Ogawa M Lawrence NA Williams LN Eggington JM Singh M Smith RA Preston BD Mutator suppression and escape from replication error-induced extinction in yeast PLoS Genet 2011 7 10 e1002282 10.1371/journal.pgen.1002282 22022273
Herr AJ, Ogawa M, Lawrence NA, Williams LN, Eggington JM, Singh M, Smith RA, Preston BD (2011) Mutator suppression and escape from replication error-induced extinction in yeast. PLoS Genet 7(10):e100228222022273 10.1371/journal.pgen.1002282
Herr AJ Kennedy SR Knowels GM Schultz EM Preston BD DNA replication error-induced extinction of diploid yeast Genetics 2014 196 3 677 691 10.1534/genetics.113.160960 24388879
Herr AJ, Kennedy SR, Knowels GM, Schultz EM, Preston BD (2014) DNA replication error-induced extinction of diploid yeast. Genetics 196(3):677–69124388879 10.1534/genetics.113.160960
Kesten H Stigum BP Limit theorems for decomposable multi-dimensional Galton-Watson processes J Math Anal Appl 1967 17 2 309 338 10.1016/0022-247X(67)90155-2
Kesten H, Stigum BP (1967) Limit theorems for decomposable multi-dimensional Galton-Watson processes. J Math Anal Appl 17(2):309–33810.1016/0022-247X(67)90155-2
Morrison A Johnson AL Johnston LH Sugino A Pathway correcting DNA replication errors in saccharomyces cerevisiae EMBO J 1993 12 4 1467 1473 10.1002/j.1460-2075.1993.tb05790.x 8385605
Morrison A, Johnson AL, Johnston LH, Sugino A (1993) Pathway correcting DNA replication errors in saccharomyces cerevisiae. EMBO J 12(4):1467–14738385605 10.1002/j.1460-2075.1993.tb05790.x
Mullikin TW Limiting distributions for critical multitype branching processes with discrete time Trans Am Math Soc 1963 106 3 469 494 10.1090/S0002-9947-1963-0144386-6
Mullikin TW (1963) Limiting distributions for critical multitype branching processes with discrete time. Trans Am Math Soc 106(3):469–49410.1090/S0002-9947-1963-0144386-6
Nicholson MD, Cheek D, Antal T (2022) Mutation accumulation in exponentially growing populations. arXiv:2208.02088
Ogura Y Asymptotic behavior of multitype Galton–Watson processes J Math Kyoto Univ 1975 15 2 251 302
Ogura Y (1975) Asymptotic behavior of multitype Galton–Watson processes. J Math Kyoto Univ 15(2):251–302
Olver FWJ, Olde Daalhuis AB, Lozier DW, Schneider BI, Boisvert RF, Clark CW, Miller BR, Saunders BV, Cohl HS, McClain MA (Eds). NIST Digital Library of Mathematical Functions. https://dlmf.nist.gov/, Release 1.1.10 of 2023-06-15
Schumacher TN Schreiber RD Neoantigens in cancer immunotherapy Science 2015 348 6230 69 74 10.1126/science.aaa4971 25838375
Schumacher TN, Schreiber RD (2015) Neoantigens in cancer immunotherapy. Science 348(6230):69–7425838375 10.1126/science.aaa4971
Schumacher TN Scheper W Kvistborg P Cancer neoantigens Annu Rev Immunol 2019 37 173 200 10.1146/annurev-immunol-042617-053402 30550719
Schumacher TN, Scheper W, Kvistborg P (2019) Cancer neoantigens. Annu Rev Immunol 37:173–20030550719 10.1146/annurev-immunol-042617-053402
Sevast’yanov Boris Alexandrovich Transient phenomena in branching stochastic processes Theory Prob Appl 1959 4 2 113 128 10.1137/1104011
Sevast’yanov Boris Alexandrovich (1959) Transient phenomena in branching stochastic processes. Theory Prob Appl 4(2):113–12810.1137/1104011
Soriano I Vazquez E De Leon N Bertrand S Heitzer E Toumazou S Bo Z Palles C Pai CC Humphrey TC Expression of the cancer-associated DNA polymerase ε p286r in fission yeast leads to translesion synthesis polymerase dependent hypermutation and defective DNA replication PLoS Genet 2021 17 7 e1009526 10.1371/journal.pgen.1009526 34228709
Soriano I, Vazquez E, De Leon N, Bertrand S, Heitzer E, Toumazou S, Bo Z, Palles C, Pai CC, Humphrey TC et al (2021) Expression of the cancer-associated DNA polymerase p286r in fission yeast leads to translesion synthesis polymerase dependent hypermutation and defective DNA replication. PLoS Genet 17(7):e100952634228709 10.1371/journal.pgen.1009526
Summers Jesse Litwin Samuel Examining the theory of error catastrophe J Virol 2006 80 1 20 26 10.1128/JVI.80.1.20-26.2006 16352527
Summers Jesse, Litwin Samuel (2006) Examining the theory of error catastrophe. J Virol 80(1):20–2616352527 10.1128/JVI.80.1.20-26.2006
Tejero H, Montero F, Nuño JC (2016) Theories of lethal mutagenesis: from error catastrophe to lethal defection. Quasispecies Theory Exp Syst, 161–179
Tilk S Tkachenko S Curtis C Petrov DA McFarland CD Most cancers carry a substantial deleterious load due to Hill–Robertson interference Elife 2022 11 e67790 10.7554/eLife.67790 36047771
Tilk S, Tkachenko S, Curtis C, Petrov DA, McFarland CD (2022) Most cancers carry a substantial deleterious load due to Hill–Robertson interference. Elife 11:e6779036047771 10.7554/eLife.67790
Tjørve KMC Tjørve E The use of Gompertz models in growth analyses, and new Gompertz-model approach: an addition to the unified-richards family PloS One 2017 12 6 e0178691 10.1371/journal.pone.0178691 28582419
Tjørve KMC, Tjørve E (2017) The use of Gompertz models in growth analyses, and new Gompertz-model approach: an addition to the unified-richards family. PloS One 12(6):e017869128582419 10.1371/journal.pone.0178691
Topatana W Juengpanich S Li S Cao J Hu J Lee J Suliyanto K Ma D Zhang B Chen M Advances in synthetic lethality for cancer therapy: cellular mechanism and clinical translation J Hematol Oncol 2020 13 1 22 10.1186/s13045-020-00956-5 31900191
Topatana W, Juengpanich S, Li S, Cao J, Hu J, Lee J, Suliyanto K, Ma D, Zhang B, Chen M et al (2020) Advances in synthetic lethality for cancer therapy: cellular mechanism and clinical translation. J Hematol Oncol 13:1–2231900191 10.1186/s13045-020-00956-5
Vignuzzi M Stone JK Andino R Ribavirin and lethal mutagenesis of poliovirus: molecular mechanisms, resistance and biological implications Virus Res 2005 107 2 173 181 10.1016/j.virusres.2004.11.007 15649563
Vignuzzi M, Stone JK, Andino R (2005) Ribavirin and lethal mutagenesis of poliovirus: molecular mechanisms, resistance and biological implications. Virus Res 107(2):173–18115649563 10.1016/j.virusres.2004.11.007
