Universal Asymptotic Clone Size Distribution for General Population Growth

Deterministically growing (wild-type) populations which seed stochastically developing mutant clones have found an expanding number of applications from microbial populations to cancer. The special case of exponential wild-type population growth, usually termed the Luria–Delbrück or Lea–Coulson model, is often assumed but seldom realistic. In this article, we generalise this model to different types of wild-type population growth, with mutants evolving as a birth–death branching process. Our focus is on the size distribution of clones—that is the number of progeny of a founder mutant—which can be mapped to the total number of mutants. Exact expressions are derived for exponential, power-law and logistic population growth. Additionally, for a large class of population growth, we prove that the long-time limit of the clone size distribution has a general two-parameter form, whose tail decays as a power-law. Considering metastases in cancer as the mutant clones, upon analysing a data-set of their size distribution, we indeed find that a power-law tail is more likely than an exponential one.


Introduction
Cancerous tumours spawning metastases, bacterial colonies developing antibiotic resistance or pathogens kickstarting the immune system are examples in which events in a primary population initiate a distinct, secondary population. Regardless of the scenario under consideration, the number of individuals in the secondary population, and how they are clustered into colonies, or clones, is of paramount importance. An approach which has offered insight has been to bundle the complexities of the initiation process into a mutation rate and assume that the primary, or wild-type, population seeding the secondary, or mutant, population is a random event.
This method was pioneered by microbiologist Salvador Luria and theoretical physicist Max Delbrück (Luria and Delbrück 1943). In their Nobel prize winning work, they considered an exponentially growing, virus susceptible, bacterial population. Upon reproduction, with small probability, a virus resistant mutant may arise and initiate a mutant clone. This model was contrasted with each wild-type individual developing resistance upon exposure to the virus with a constant probability per individual. By considering the variance in the total number of mutants in each case, they demonstrated that bacterial evolution developed spontaneously as opposed to adaptively in response to the environment.
In the original model of Luria and Delbrück, both wild-type and mutant populations grow deterministically, with mutant initiation events being the sole source of randomness. Lea and Coulson (1949) generalised this process by introducing stochastic mutant growth in the form of the pure birth process and were able to derive the distribution of the number of mutants for neutral mutations. This was again extended by Bartlett (1955) and later Kendall (1960), who considered both populations developing according to a birth process. An accessible review discussing these formulations is given by Zheng (1999).
Recent developments have focused on cancer modelling, where usually mutant cell death is included in the models. The main quantity of interest in these studies has been the total number of mutant cells. Explicit and approximate solutions appeared for deterministic, exponential wild-type growth, corresponding to a fixed size wild-type population (Angerer 2001;Dewanji et al. 2005;Iwasa et al. 2006;Komarova et al. 2007;Keller and Antal 2015), and fully stochastic wild-type growth either at fixed time or fixed size (Durrett and Moseley 2010;Antal and Krapivsky 2011;Kessler and Levine 2015). An exciting recent application has been to model emergence of resistance to cancer treatments (Kessler et al. 2014;Bozic et al. 2013;Bozic and Nowak 2014). The current study continues in this vein with our inspiration being primary tumours (wild-type) seeding metastases (mutant clones).
Interestingly, in the large time small mutation rate limit, the clone size distribution at a fixed wild-type population size coincides for stochastic and deterministic exponential wild-type growth (Kessler and Levine 2015;Keller and Antal 2015). The intuition behind this observation is that a supercritical birth-death branching process converges to exponential growth in the large time limit, and, for a small mutation rate, mutant clones are initiated at large times. So asymptotically the two methods are equivalent, but the deterministic description of the wild-type population has twofold advantages: (i) the calculations are much simpler in this case (Keller and Antal 2015), and (ii) the method can be easily generalised to arbitrary growth functions. This is the programme that we develop in the present paper.
The present work differs from previous approaches in two ways. Firstly, motivated by populations with environmental restrictions, we move away from the assumption of exponential wild-type growth, a setting which has received limited previous consideration as discussed in Foo and Michor (2014). We shall first review and extend results for the exponential case and then provide explicit solutions for power-law and logistic growth. Next, we present some general results which are valid for a large class of growth functions. This extends the classic results found in Kendall (1948), Athreya and Ney (2004), Karlin and Taylor (1981), Tavare (1987) and recent work in Tomasetti (2012), Houchmandzadeh (2015) who considered the wild-type population growth rate to be time-dependent but coupled with the mutant growth rate. Secondly, rather than the total number of mutants, our primary interest is on the distribution of mutant number in the clones initiated by mutation events. This complements Hanin et al. (2006), which allowed deterministic wild-type and mutant growth, and the treatment of clone sizes for constant wild-type populations found in Dewanji et al. (2011). While we focus on clone sizes, we demonstrate that the distribution for the total number of mutants follows as a consequence, and hence, results hold in that setting also.
The outline of this work is as follows. We define our model in Sect. 2, utilising formalism introduced in Karlin and Taylor (1981), and demonstrate a mapping between the mutant clone size distribution and the distribution for the total number of mutants. The exact time-dependent size distribution is given for exponential, power-law and logistic wild-type growth function in Sect. 4. Section 5 pertains to universal features of the clone size distribution and contains our most significant results. There, for a large class of wild-type growth functions, we demonstrate a general two-parameter distribution for clone sizes at large times. The distribution has power-law tail behaviour which corroborates previous work (Iwasa et al. 2006;Durrett and Moseley 2010;Williams et al. 2016). Large time results are also given for the mean and variance of the clone sizes under general wild-type growth. Adopting the interpretation of the wild-type population as the primary tumour and mutant clones as metastases, we test our results regarding the tail of the distribution on empirical metastatic data in Sect. 6. Section 7 considers alternative methods to ours, and we give some concluding remarks in Sect. 8.

Model
In our model, a wild-type population gives rise to mutants during reproduction events. The arisen mutant also reproduces, and so mutant clones stem from the original initiating mutant's progeny. In many applications, the wild-type population is significantly larger than the mutant clones, and so we treat the wild-type population's growth as deterministic, with size dictated by a time-dependent function n τ . The mutant clones are smaller in comparison, and so their growth is stochastic. For logistic wild-type growth, a sample realisation of the process is shown in Fig. 1. The exact formulation is now given. for deterministic logistic wild-type growth, with a carrying capacity of 50, and stochastic mutant growth. Note that we typically assume the wild-type population is much larger than individual clones

The Birth-Death Process
Stochastic growth of mutants will follow a birth-death branching process (Athreya and Ney 2004). Time is scaled such that each mutant has unit birth rate and death rate β. A brief note on converting our results to the case when the birth rate is arbitrary is given in "Appendix B". Let Z t be the size of a population at time t, with Z 0 = 1. The forward Kolmogorov equation for the distribution is given by (1) with k ≥ 1. Its solution in terms of the generating function, given on page 76 of Bartlett (1955), is Due to our timescale, β is the probability of eventual extinction for a mutant clone for β ≤ 1, and λ is the mutant fitness. When β = 0, and so the stochastic proliferation follows a pure birth or Yule process, the mutants will be denoted immortal. By expanding the generating function around s = 0, we obtain for the probability of the population size being k a geometric distribution with a modified zero term with the shorthand notation For the particular case of a critical branching process, i.e. when birth and death rates are equal, the above probabilities are simplified by observing

Mutant Clone Size Distribution
Here, we employ standard methods as outlined in, for instance, Karlin and Taylor (1981), Dewanji et al. (2005). The system is observed at a fixed time t, and we let the number of wild-type individuals be denoted by n τ for 0 ≤ τ ≤ t. Since mutants are produced by wild-type individuals, the rate of mutant clone initiations will be proportional to the product of n τ and the mutation rate μ. More precisely, the process of clone initiations is an inhomogeneous Poisson process (Karlin and Taylor 1998) with intensity μn τ . Let the Poisson random variable K t denotes the number of clones that have been initiated by t, which has mean Now, assuming K t > 0, we consider a mutant clone sampled uniformly from the K t initiated clones and denote its size to be the random variable Y t . The clone was initiated at the random time T , and as we must have T ≤ t, the density of T is given by Where is the expected number of clones seeded when the mutation rate is unity. The size of the clone is dictated not only by the initiation time but also by its manner of growth, here the birth-death process. Hence, by conditioning on the arrival time, we have An immediate consequence is that the generating function of the clone size is given by where Z t (s) is the generating function of the birth-death process (2).
We make the following remarks on the above. (i) The mutation rate μ does not appear in the density for initiation times in (6); hence mutant clone sizes are independent of the mutation rate and thus all following results regarding clone sizes will be also. (ii) The integral in (8) is a convolution, and as convolutions commute, we may swap the arguments of the integrand functions (n τ Z t−τ ↔ n t−τ Z τ ). (iii) If we start with n 0 wild-type individuals, so the wild-type follows m τ = n 0 n τ , then both the numerator and denominator in (6) will have a factor of n 0 , which cancel. So henceforth, apart from when n 0 = 0 (used occasionally for analytic convenience), we set n 0 = 1 without loss of generality. (iv) By similar logic, a positive random amplitude for the wild-type growth function, i.e. m τ = Xn τ for a general positive random variable X , would also cancel, and so our results on clone sizes hold in that case also.

Mapping Distributions: Clone Size to Total Mutant Number
This section is related to the classic Luria-Delbrück problem. Let B t be the total number of mutants existing at time t. Then, B t is the sum of K t generic clones where all (Y t ) i are iid random variables specifying the clone sizes. As such, B t is a compound Poisson random variable, and hence its generating function is which can be derived by conditioning on K t . It follows that The link between the mass functions of the mutant clone size, Y t , and the total number of mutants, B t , is given by the recursion This relationship may be found as Lemma 2 in Zheng (1999), and a short proof is provided for convenience in "Appendix B", Lemma B1. Hence, while we may initially work in the setting of size distribution of a single clone, by the above discussion, results are transferable to the total number of mutants case. Often long-time results are sought, which significantly reduces the complexity of the distributions. For any fixed positive mutation rate, in the long-time limit, an infinite number of clones will have been initiated, and thus, the probability distributions of B t will not be tight (Durrett 1996). A common solution to this problem is the Large Population-Small Mutation limit (Keller and Antal 2015), where θ = μn t is kept constant. Then, for exponential wild-type growth, n τ = e δτ , (or exponential-type, see Sect. 5), the expected number of initiated clones, E(K t ), tends to θ/δ for large times. Hence, we see that demonstrating that the limit of the clone size distribution is of primary concern. Furthermore, if the expected number of initiated clones is small, we have the following proposition, whose proof can be found in "Appendix B".
Proposition 1 For a small expected number of initiated clones, conditioned on survival, the size of a single clone and the total number of mutants are approximately equal in distribution. That is, One immediate consequence of this result is that for immortal mutants (β = 0) and This agrees with intuition as for small enough E(K t ), we expect only 0 or 1 clones to be initiated, and hence, the total number of mutants will be dictated by the clone size distribution. With exponential wild-type growth, this approximation was used in Iwasa et al. (2006) to investigate drug resistance in cancer.

Finite Time Clone Size Distributions
Three particular cases of wild-type growth function, n τ , will be considered in detail, namely: exponential, power-law and logistic (Fig. 2). Exponential and logistic growth are widely used in biological modelling (Murray 2002). For the power-law cases, under the assumption that the radius of a spherical wild-type population is proportional to time, quadratic and cubic power-law growth represents mutation rates proportional to the surface area and volume, respectively. In each case, we give the generating function and probability mass function. We stress again that the mutation rate and an arbitrary positive prefactor for n τ cancel in (8) and so are irrelevant for our results.

Exponential Wild-Type Growth
Let the wild-type population grow exponentially, that is n τ = e δτ with δ > 0 and so from (7), a t = e δt −1 δ . The distribution for the total number of mutants, B t , was  (2015), and we follow their notation by letting γ = δ/λ. Using (10) and the results found in section 3 of Keller and Antal (2015), the generating function is Similarly, the mass function is Here, F a, b c ; z is Gauss's hypergeometric function, and (a) k is the Pochhammer symbol defined in "Appendix A". The above expressions are given in terms of n t to allow easy comparison to the formulas in Keller and Antal (2015). For these exact timedependent formulas, their form is somewhat cumbersome; however, simpler long-time limit expressions are given in Sect. 5. A reduction in complexity is also obtained for the special case of neutral mutants (δ = λ) where, by using (24), the generating function in (12) simplifies to If additionally the neutral mutants are immortal, the above expression further simplifies to The probabilities are then concisely given by which corresponds to the classical Lea-Coulson result (Lea and Coulson 1949)

Power-Law Wild-Type Growth
Now, we assume that the wild-type population grows according to a general power-law n τ = τ ρ , for some non-negative integer ρ, and therefore, a t = t ρ+1 ρ+1 . With power-law wild-type growth and stochastic mutant proliferation, the mutant clone size generating function is given by Here, Li i (s) is the polylogarithm of order i defined in "Appendix A". Details of the derivation are given in "Appendix C". For immortal mutants, the mass function may be explicitly written as If mutants may die, the exact mass function is most easily obtained via Cauchy's integral formula which may be efficiently computed using the fast Fourier transform. For a brief discussion on implementation, see Antal and Krapivsky (2010) and references therein.
Note for ρ ≥ 1, n 0 = 0 which, while useful for analytic tractability, is unrealistic. This can be overcome by letting n τ = n 0 + τ ρ . Then, by splitting the integral in the generating function (9) and using the above analysis, one can obtain the mass function for any n 0 . However, for practical purposes, the contribution of n 0 is negligible.

Constant Size Wild-Type
For the specific power-law growth when ρ = 0, i.e. n τ = 1 (recall that this is equal to the general case when n τ = n 0 ), we recover some classical results for constant immigration (Kendall 1948). We note that the distribution of the ordered clone size, depending on initiation time, was discussed in Jeon et al. (2008). From (13) with ρ = 0, the generating function is with S t as given in (4). By expanding this generating function in terms of s we obtain the probabilities Then, using (10) with the clone sizes (15) we obtain the generating function of the total number of mutants and from the binomial theorem we also get the probabilities We recognise this as a negative binomial distribution under the interpretation that B t is the number of failures until μ successes, with failure probability S −1 t . This result for B t was first derived by Kendall (1948) who was attempting to explain the appearance of the logarithmic distribution for species number when randomly sampling heterogeneous populations, conjectured by R.A. Fisher. From the distribution of B t , by an argument which may be considered a special case of Proposition 1, he derived that for constant rate initiation, the clone size conditioned on non-extinction is logarithmically distributed again with parameter S −1 t , which can be obtained via (16). Constant immigration may imply a constant size source; hence, mutants with equal birth and death rates (i.e. evolving as a critical branching process) are particularly interesting. This case yields analogous formulas to those above but S t is replaced with the expression given in (5).

Logistic Wild-Type Growth
Starting from a population of one and having a carrying capacity K , logistic growth is given by n τ = K e λτ K +e λτ −1 . We assume neutral mutations, i.e. λ is also the wild-type growth rate. Integrating the growth function gives a t = K λ log e λt n t . We aim to calculate the generating function using (9). Recalling the definition of Z t−τ (s) we observe that where C is an integration constant. Therefore, the generating function is .
Agreeing with intuition for K = 1, we recover the generating function of the constant case, and lim K →∞ Y t (s) gives the generating function for exponential wild-type growth. Therefore, the logistic case interpolates between the constant and exponential growth cases. The mass function can be obtained by expanding the non-logarithmic and logarithmic function in Y t (s) and using the Cauchy product formula. However, this method provides little insight, and numerically, it is simpler to use the fast Fourier transform.

Monotone Distribution and Finite Time Cut-Off
We conclude this section by demonstrating general features that exist in the clone size distribution at finite times. Again proofs are provided in "Appendix C". Firstly, we see that, regardless of the particular wild-type growth function, the monotone decreasing nature of the mass function for the birth-death process is preserved in the clone size distribution.

Proposition 2 As long as n
Whether P(Y t = 0) ≥ P(Y t = 1) depends on n τ and t, but the inequality is typically true for long times. Note that in contrast, the mass function of the total number of mutants is not monotone in general (Keller and Antal 2015). Now restricting ourselves to the λ > 0 case, as an example, consider the mass function when the size of the wild-type population is constant, which is given by (16), and specifically for k ≥ 1. For any moderate t, S −1 t is typically close to unity but for large k, S −k t will become the dominant term in the mass function, dictating exponential decay. We term this a cut-off in the distribution which occurs at approximately k = O(e λt ). It is an artefact of the mass function for the birth-death process (3). Hence, we will have (at least) two behaviour regimes for the mass function for finite times. Here, we show that this cut-off exists generally for finite times.
Note that S t > 1 for finite t, and S t → 1 exponentially fast for large t. Hence, the cut-off will disappear for long times and the subexponential factor, discussed in detail in Sect. 5, will completely determine the tail of the distribution. Also notice that the power-law cases, n τ = τ ρ , for ρ ≥ 1 are not covered as, to make the analysis tractable, they artificially start at n 0 = 0. However, the generating function in this case (13) also has its closest to origin singularity at S t so the cut-off exists there also.

Universal Large Time Features
Here, we give results regarding the large time behaviour of our model which is relevant in many applications and also provides general insight. In many applications, the cutoff location (k = O(e λt )) is so large that the distribution at or above this point is of little relevance, and hence, for this purpose the limiting approximations now discussed are of particular interest. Using the notation of Theorem 1, this section investigates the large time form of Θ t (k). The proofs for the results presented in this section can be found in "Appendix D". We highlight the power-law decaying, "fat" tail found in each case. Henceforth, we again assume λ > 0, i.e. a supercritical birth-death process.

General Wild-Type Growth Functions
To give general results, we introduce the following assumption which will be assumed to hold for the remainder of this section.
Assumption 1 For wild-type growth function n τ , we assume (i) n τ = 0 for τ < 0, continuous for τ > 0 and right continuous at τ = 0. (ii) n τ is positive and monotone increasing for τ > 0. (iii) For x ≥ 0 the limit lim t→∞ n t−x /n t exists, is positive and finite.
We note that the cases discussed in Sect. 4 are all covered by Assumption 1. The reason for the seemingly arbitrary limit assumed in (iii) becomes clear with the following result which is an application of the theory of regular variation found in Bingham et al. (1987). Often the long-time behaviour of the clone size distribution may be separated into δ * > 0 and δ * = 0, and so we give the following definition (Flajolet and Sedgewick 2009).

Definition 1 Consider a real valued function
Simple examples of subexponential functions are e √ t , log(t), t ρ , while e δt , e δt t ρ are of exponential-type, with δ, ρ ∈ R.

Mean and Variance
We now address the asymptotic properties of the clone size distribution by first discussing its mean and variance.

Theorem 2 With s i (t) subexponential functions such that s
as t → ∞.

Large Time Clone Size Distribution
Turning to the distribution function, we have the following result regarding the generating function at large times.
This result is made clearer in the next corollary, in which the cases of exponentialtype and subexponential growth are separated. This is as, for δ * > 0, lim t→∞ n t a t → δ * .
For a proof, see Lemma D1. Consequently, in the exponential-type setting, the limiting result is a proper probability distribution, while in the subexponential case it is not. We can interpret this as the clone sizes staying finite in the exponential case but grow to infinity for subexponential cases at large times. Henceforth, for brevity, we do not impose such a separation but the reader should note that for exponential-type growth the above limit holds and may simplify further results.
where the second expression is the γ * → 0 limit of the first expression. Then for t → ∞ the probabilities for exponential-type growth γ * > 0 are and for subexponential growth (γ * = 0) This result is exemplified in Fig. 4. The expressions obtained in the δ * > 0 case also appeared as an approximation in Kessler and Levine (2015) for the total number of mutants with stochastic wild-type and mutant growth when the mean number of clones is small. This can now be interpreted as an application of Proposition 1.
The case of immortal mutants does not simplify the above expressions for subexponential growth, but for exponential-type growth, by applying (23) then (22) to the limiting generating function, we have the following link to the Yule-Simon distribution which appears often in random networks (Simon 1955;Krapivsky and Redner 2001). Transition to the asymptotic regime as described in Corollary 1. For subexponential wild-type growth, the mass functions tend to k −1 behaviour, while for exponential-type it tends to k −1−γ * . Here, t = 20 and all other parameters are as given in Fig. 2 Corollary 2 For immortal mutants with exponential-type wild-type growth the clone size distribution Y t follows a Yule-Simon distribution with parameter δ * for large times. That is, for With immortal, neutral (δ * = 1) mutants we have which is in agreement with the long-time limit of (4.1). For immortal mutants and exponential-type growth, as the clone size distribution tends to a Yule-Simon distribution, we expect power-law tail behaviour at large times (Newman 2005). Interestingly, we see that this behaviour holds when we include mutant death and have general wild-type growth.
Corollary 3 At large times, the tail of the clone size distribution follows a power-law with index 1 + γ * . More precisely,

Large Time Distribution for Total Number of Mutants
Finally, to conclude this section, we give the corresponding results for the total number of mutants B t in the often used Large Population-Small Mutation limit.

Theorem 4 Letting θ = μn t be constant and with s i (t) subexponential functions as in Theorem
and we have the following tail result

Tail Behaviour in Empirical Metastatic Data
Given the above discussion we expect, for a large class of wild-type growth functions, to see power tail behaviour on approach to the exponential cut-off in the clone size distribution. We take the first steps to verify this theoretical hypothesis by analysing an empirical metastatic data. In this setting, the wild-type population is the primary tumour and mutant clones are the metastases. Our data are sourced from the supplementary materials in Bozic et al. (2013). These data are taken from 22 patients; 7 with pancreatic ductal adenocarcinomas, 11 with colorectal carcinomas, and 6 with melanomas. One patient had only a single metastasis so we discard this data. Of the 21 remaining patients, the number of cells in a single metastasis ranged from 6 × 10 6 to 2.23 × 10 9 . Our theoretical model predicts a cut-off in the distribution around k = e λt . Taking some sample parameters from the literature, namely λ = 0.069/day (Diaz et al. 2012), and t = 14.1 years (Yachida et al. 2010), this leads to a cut-off around k ≈ 10 154 cells. Due to the enormity of this value, we ignore the cut-off here. Additionally, as the minimum observed metastasis size is 6 × 10 6 cells, we assume that all data points are sampled from the tail of the distribution.
For each of the data-sets, we examine the likelihood ratio to determine whether the data is more likely sampled from a power-law decaying or geometrically decaying distribution. Nineteen of the 21 data-set return the power-law hypothesis as more plausible which is in agreement with the theoretical prediction. Both are single parameter distributions, and maximum likelihood analysis was utilised to estimate the parameters. The methodology outlined in Clauset et al. (2009) was broadly followed, and Hence, values to left of figures are more significant. a Likelihood ratiô R for each data-set. Points above the horizontal line suggest data-set is from a power-law distribution over a geometric distribution. b Estimateω for each data-set, determined via maximum likelihood. c Normalised log-likelihood function for best data-set. Vertical bars show the likelihood interval brief details regarding calculating maximum likelihood estimates (MLEs) are given in "Appendix E". We note that in this context the likelihood ratio point esimator returns equivalent results to the Akaike information criterion widely used in model selection (Burnham and Anderson 1998). Under the power-law model, P(Y t = k) ∝ k −ω , for 20 of the 21 data-sets, we find the point estimate of the power-law index,ω, lies in [−2, −1]. The outlier comes from the smallest data-set (3 metastases). Due to the small size of data-sets, we recognise the influence of statistical fluctuations. Details of the likelihood ratio are as follows. Let y = (y 1 , . . . , y N ) be a dataset of size N . We test the hypothesis that y is drawn from a power-law distribution, P 1 (Y t = k) = C 1 k −ω , versus that it is sampled from a geometric distribution, P 2 (Y t = k) = C 2 p(1 − p) k , where C 1 , C 2 are normalising constants and p is the parameter for the geometric distribution. The log-likelihood ratio iŝ whereR > 0 gives support to the hypothesis that the data is drawn from the powerlaw distribution with MLE exponentω, over the geometric distribution with MLE parameterp. The results are given in Fig. 5a.
Assuming a power-law distribution, the maximum likelihood estimates for the exponent ω for each data-set are given in Fig. 5b. Due to the small sample size of our data-sets and the high variance in the distribution, we do not derive confidence intervals via normal distribution approximations. Instead we show the normalised log-likelihood, log L(ω)/L(ω), for our best data-set, with N = 30, in Fig. 5c, where L(ω) is the likelihood function. Also, following Hudson (1971), we demonstrate the likelihood interval defined as If a large sample size was possible this interval would correspond to a 95.4% confidence interval. For the data-set with N = 30, we numerically determined I (ω) = [1.295, 1.616], demonstrated as the domain between the vertical bars in Fig. 5c.

Deterministic Approximation
In order to circumvent the complexity introduced by the birth-death process, one might be tempted to simply assume the mutant clone size grows according to e λτ , the mean of the birth-death process. This approach corroborates our results regarding the tail of the size distribution. Indeed, the clone size density may be found to be which has support [1, e λt ]. This formula can also be found in Hanin et al. (2006). Then, as in Sect. 5 under Assumption 1, Thus, asymptotically the density has the same behaviour as the tail of the limiting result given in Corollary 3, but with a different amplitude. However, despite this agreement, the densities given by (17) for specific wild-type growth function differ significantly compared with stochastic mutant proliferation. Letting Y Stoch t be the clone size distribution with stochastic mutant growth and Y Det t be its deterministic approximation specified by (17), we may quantify the approximation error, at least for the moments, by the following theorem, whose proof can be found in "Appendix F".

Time-Dependent Rate Parameters
Some authors Houchmandzadeh (2015), Tomasetti (2012) have previously considered the case where all rates in the system are multiplied by a time-dependent function, say z(τ ). This is relevant in the scenario where both the wild-type and mutant populations have their growth restricted simultaneously by environmental factors, for example exposure to a chemotherapeutic agent. We observe that under a change of timescale this system is equivalent to our setting with exponential wild-type growth. This is due to the following argument.
In this setting, the wild-type population is governed by Mutant clones are now initiated at a rate μz(τ )n τ . Let Z t be the size of a mutant population governed by the birth-death process with time-dependent rates. Once initiated, the size distribution obeys the forward Kolmogorov equation for time-dependent stochastic mutant proliferation If we let then under a new timescale, τ = F −1 (τ ) , the mutant clone initiation will occur at a rate μn τ . Further, using the chain rule to express (18) and (19) in terms of τ , we see that n τ = e λτ and that the forward Kolmogorov equation (19) becomes (1). Thus, under a time-rescaling, all dynamics are equivalent to the system with exponential wild-type growth and stochastic mutant proliferation with constant birth and death rates, as studied in this article or in Keller and Antal (2015).

Poisson Process Characterisation of Tail
Complementing Corollary 3 in Sect. 5, following Tavare (1987), we can also describe the size distribution for large clones at long times via a Poisson process in the following way. Let (Z (i) (t)) i≥1 be independent copies of the birth-death process as in Sect. 2 and (T i ) i≥1 ⊂ (0, ∞) be the points of a of Poisson process with intensity μn τ , for τ ≥ 0. The T i represent the clone arrival times, and so K t is the number of T i less than or equal to t. Let us consider the size of the first clone. By known results about the large time behaviour of the birth-death process (Athreya and Ney 2004), as t → ∞, The distribution of the limiting random variable W 1 is composed of a point mass at 0 and an exponential random variable, precisely Analogously, with the details given in Tavare (1987) (Theorem 3), the limiting behaviour of the time-ordered clone sizes is given by where W 1 is as before and all W i are iid. The random sequence (e −λT i W i ) i≥1 takes nonnegative real values; however, if we restrict our attention to only the positive elements (that is clones that do not die), then these can be taken to be points from a nonhomogeneous Poisson process. More precisely, the set {σ j } j≥1 := {e −λT i W i } i≥1 \ {0} are the points (in some order) from a Poisson process on (0, ∞) with mean measure The proof of the above only requires minor modification from that of Theorem 4 in Tavare (1987). The Poisson process description of the large clones, at large times, can also offer insight into further properties of the system, including links to the Poisson-Dirichlet distribution, see Tavare (1987), Durrett (2015). With regards to the present article, the interesting point is that for fixed ε > 0, as the number of σ j > ε is finite almost surely, we may sample unformly from this set (i.e. {σ j } j≥1 ∩ (ε, ∞)) and construct a random variable Y ε with distribution (20). The new variable Y ε can be related to the previously considered random variable Y t by the following result, whose proof is contained in "Appendix F".
Theorem 6 For ε > 0, with Y ε as above, Of note is the reappearance of power-law behaviour with a cut-off in the density of Y ε . For example in the constant wild-type case, n τ = 1, the density, using (20), is given by For exponential growth with neutral mutants, n τ = e λτ , Note that the exponents in the power-law terms is equal to that given in Corollary 3, indicating the two approaches are complimentary.
tion is dictated by a generic deterministic growth function, and the mutant growth is stochastic. This shifts the focus from previous studies which have mostly been concerned with exponential, or mean exponential, wild-type growth, and considered the total number of mutants. Results for the total number of mutants are, however, easily obtained from the clone size distribution. The special cases of exponential, power-law and logistic wild-type growth were treated in detail, due to their extensive use in models for various applications. Utilising a generating function centred approach, exact time-dependent formulas were ascertained for the probability distributions in each case. Regardless of the growth function, the mass function is monotone decreasing and the distribution has a cut-off for any finite time. This cut-off goes to infinity for large times and is often enormous in practical applications; hence, we focused on the approach to the cut-off.
We found that the clone size distribution behaves quite distinctly for exponentialtype versus subexponential wild-type growth. Although the probability of finding a clone of any given size stays finite as t → ∞ for exponential-type growth, it tends to zero for subexponential type. Despite these differences, with a proper scaling, for a large class of growth functions, we proved that the clone size distribution has a universal long-time form. This long-time form possesses a power-law "fat" tail which decays as 1/k for subexponential wild-type growth, but faster for exponential-type growth. This can be intuitively understood as the tail distribution represents clones that arrive early, and the chance that a clone is initiated early in the process is larger for a slower growing wild-type function. Hence, we expect a "fatter" tail in the subexponential case.
Note that although we consider the case of subexponential wild-type growth, surviving mutant clones will grow exponentially for large time, which can be unrealistic in some situations. Stochastic growth which accounts for environmental restrictions, for instance the logistic branching process, introduces further technical difficulties and is left for future work. We do note that, despite the drawbacks of deterministic mutant growth as discussed in Sect. 7, when both the wild-type and mutant populations grow deterministically as τ ρ , it is easy to see that for large times the clone size distribution still displays a power-law tail, lim t→∞ t f Y t (y) = ρ+1 ρ y 1/ρ−1 . An underlying motivation for this work is the scenario of primary tumours spawning metastases in cancer. We test our hypothesis regarding a power-law tail in metastasis size distributions by analysing empirical data. For 19 of 21 data-sets, the power-law distribution is deemed more likely than an exponentially decaying distribution. The exponent of the power-law decay was estimated in each case and found to lie between −1 and −2. Interpreting this in light of our theory, either the primary tumour had entered a subexponential growth phase or, if one assumes exponential primary growth, the metastatic cells had a fitness advantage compared to those in the primary. Either way we can conclude that, for the majority of patients, the metastases grew faster than the primary tumour.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.
For any analytic function f (z) = n≥0 f n z n , we denote the nth coefficient as Theorem A1 (Flajolet and Sedgewick 2009: Exponential Growth Formula) If f (z) is analytic at 0 and R is the modulus of a singularity nearest the origin in the sense that R := sup{r ≥ 0| f is analytic in |z| < r }. Then the coefficient We utilise several results from Bingham et al. (1987) on the theory of regularly varying functions which we now define.
Definition 2 (Bingham et al. 1987) A Lebesgue measurable function f : R + → R that is eventually positive is regularly varying (at infinity) if for some κ ∈ R, The notation f ∈ RV κ will be used and we will denote f ∈ RV 0 as slowly varying functions.
Theorem A2 (Bingham et al. 1987: Characterisation Theorem) Suppose f : R + → R is measurable, eventually positive, and lim t→∞ f (t) exists, and is finite and positive for all x in a set of positive Lebesgue measure. Then, for some κ ∈ R,

Appendix B: Proofs for Section 3
In this work, we have fixed the birth rate to be one. Other works, for example Keller and Antal (2015), use a birth-death process with birth rate α and death rate β under timescale t . Then, the timescale used in the present work is defined by t = α t . This in turn implies that all rates under t are given by dividing the corresponding rate under t by α , e.g. β = β α .

Lemma B1
Consider generating functions F(s) = n≥0 p n s n and G(s) = n≥0 q n s n where F(s) = e G(s) . Then p 0 = e q 0 and for n ≥ 1 the following recursion holds Proof Clearly, p 0 = e q 0 from F(0) = e G(0) . By differentiating F(s), we obtain F (s) = F(s)G (s), and in general which can be shown by induction using Pascal's formula for binomial coefficients.
Evaluating the above equation at s = 0 and using that F (m) (0) = m! p n and G (m) (0) = m!q n we arrive at the announced recursion.

Appendix C: Proofs for Section 4
We derive the generating function for the clone size distribution for stochastic growth and power-law wild-type growth, n τ = τ ρ , given in (13). From (9), we have It is enough to show where C is a constant of integration. This may be derived by a binomial expansion of the denominator and an identity for the incomplete gamma function, but for succinctness we simply differentiate both sides with respect to τ . First, we note that Now differentiating the right hand side of (26) yields where the equality follows by the telescoping nature of the sums. Noting that (1 − ξ e −λ(t−τ ) ) −1 = Li 0 ξ e −λ(t−τ ) + 1 and applying the limits of the integral gives the desired result.
To determine the mass function, we seek a power series representation of the generating function. We focus on the β = 0 case and thus ξ = s s−1 . By the definition of the polylogarithm and the binomial theorem Reindexing the sum, we obtain Proof of Proposition 2 Using (8), we see that for k ≥ 1 (3), it is clear that the integrand is negative for finite, positive τ giving the result.

Proof of Theorem 1
The result is an application of Theorem A1. We seek the closest to the origin singularity of which is claimed to be at S t . Indeed, we note that for |s| < S t , Z τ (s) is analytic for all τ , and as n τ is continuous, we can conclude that the I t (s) is analytic in this region also [Chapter 2, Theorem 5.4 in Stein and Shakarchi (2003)]. As n τ > 0 there exists ε > 0 such that .
The rightmost expression can be seen to have closest to origin singularity at S t and as a t Y t (s) = I t (s), by Theorem A1, we can conclude Theorem 1.

Appendix D: Proofs for Section 5
Proof of Lemma 1 Choose x ≥ 0 and let y = e t , c = e −x . Consider the function g(z) = n log(z) . Then Theorem A2(i) yields lim t→∞ n t−x n t = lim y→∞ g(yc) g(y) = c δ * = e −xδ * .
Further, Proposition A1 gives The non-negativity of δ * is dictated by the monotone increasing nature of n τ .
To prove Theorem 2, we require the following: Lemma D1 Let s 1 (t), s 2 (t) be subexponential functions, then We highlight that neither subexponential function depend on η.
Proof (i) For y = e t , n log y = g(y). Now g ∈ RV δ * hence g(y) = y δ * l(y) with l ∈ RV 0 by Theorem A2(ii). Setting s 1 (log y) = l(y), by Lemma 1, s 1 (t) is subexponential. (ii) Let g(y) be as above. With δ * > η ≥ 0 and using the change of variables s = log τ we have where the asymptotic equivalence is due to Theorem A3 applied to g(y)y −1−η ∈ RV δ * −η−1 and the final equality is by part (i). For δ * = η, by Theorem A2(ii) the integrand will be a subexponential function. Applying the same change of variables as in (27), we see by Proposition A2(i) that the integral is a slowly varying function in y and hence is subexponential in t, which we denote s 2 (t). Now for δ * < η, by Lemma 1, we may choose t large enough such that n 1/t t < e (δ * +η)/2 which by a basic result for Laplace transforms, see, e.g. Theorem 1.11 in Schiff (1999), ensures convergence to a finite, positive constant.
As an example, which will be useful for the next proof, we apply the above lemma to a t . With s 1 (t), s 2 (t) subexponential functions Proof of Theorem 2 We require the first and second moments of Z t which may be found by differentiating (2), or see Lemma F1. Then Throughout let s t be a generic subexponential function and it will be helpful to observe that the reciprocal or constant multiples of a subexponential function are subexponential. We first consider the mean. For the cases δ * = λ, applying Lemma D1(ii) to (28) with η = λ for the numerator and η = 0 for the denominator proves the claim. For δ * = λ, using Lemma D1(i) then (ii), we have That s 1 (t) diverges can be seen by applying the standard change of variables t = log(y), τ = log(s) coupled with Proposition A2(ii). Turning to the variance, with , when δ * > 2λ we apply Lemma D1(ii) term by term to (29). For δ * < λ all integrals converge and so, with C 1 , C 2 constants, The last relation is due to the monotonicity of a t and the desired representation is obtained by applying Lemma D1(ii) to a t and absorbing C 1 into s 4 (t). When λ ≤ δ * < 2λ, the same argument holds as long as we note that By Proposition A2(i) the rightmost integral is a subexponential function and as we may always choose t sufficiently large such that s Applying Lemma D1 to a −1 t t 0 n τ e −λt dτ demonstrates the contribution from the mean squared is negligible. For δ * = 2λ, we apply the same argument as in (30) to each term and this completes the proof.
In order to prove Theorem 3, we require the following lemma.

Proof of Theorem 3
To avoid division by 0 let t > 0. Taking the generating function for Y t from equation (9), we apply the change of variables u = e −λτ which gives Noting that by monotonicity n t+ log u λ n t ≤ 1, which coupled with Lemma D2 shows the integrand may be dominated. By assumption, the integrand converges, and therefore, using Lemma 1 and the dominated convergence theorem, we have lim t→∞ a t n t (Y t (s) − β) = − 1 0 u δ * /λ ξ 1 − ξ u du = −1 ξ γ * B ξ (γ * + 1, 0) The final equality follows from applying (25) then (23).

Proof of Corollary 1
The first statement is given by applying Lemma D1(ii) to a t /n t . Then, taking the limiting generating function in Theorem 3, we firstly apply (23) then (24) which yields generating function representation for γ * = 0. The mass function for γ * = 0 is simply a logarithmic expansion. For γ * > 0, we use the expression given in Appendix A of Keller and Antal (2015) to obtain a series expansion for the generating function in terms of s, and the coefficients of the expansion give the mass function.

Proof of Corollary 3
The analysis involves expanding the limiting generating function in Theorem 3 around its singularity at s = 1 and exactly mirrors that given in section 6 of Keller and Antal (2015) and so is omitted.
Proof of Theorem 4 The mean and variance can be obtained by using (11) with Theorem 2 (the second moment dominates the mean squared in all divergent cases). For the generating function, (10) gives which coupled with Theorem 3 yields the result. The map in (10) is analytic so the tail can be obtained by its expansion coupled with Corollary 3. ω = arg max ω log(L(ω)) = arg max ω log(P(Y t = y; ω)) where L(ω) is the likelihood function and the joint distribution becomes a product by independence. We derive the MLE under the assumption that the data is sampled from a distribution whose tail follows the geometric distribution, the power-law case is analogous but for the final step. Assume for at least y ≥ m, that P 2 (Y t = y; p) = C 2 p(1− p) y . Let A = P(Y t < m) (no indices are required as this quantity is assumed independent of the tail) then The log-likelihood is now given by Setting ∂ p log L( p) = 0, we solve to find the MLÊ .
The power-law case is analogous. There C 1 = 1−A ζ (ω,m) , where ζ is the Hurwitz zeta function. No closed form expression is found for the MLE and instead we have the approximation (Clauset et al. 2009) This was used for the estimates given in Sect. 6.

Appendix F: Proofs for Section 7
In order to prove Theorem 5, we require the following lemma regarding the moments of the birth-death process: Lemma Proof of Lemma F1 Recall the generating function for the birth-death process (2) whose power series representation has coefficients given by (3). The moment generating function is M Z t (s) = Z t (e s ) and hence Thus for m ≥ 1 Hence, using the second statement in Lemma F1, we have which is the desired result.
It is enough to examine the numerator. As, from (3),