Mutant Number Laws and Inﬁnite Divisibility

: Concepts of inﬁnitely divisible distributions are reviewed and applied to mutant number distributions derived from the Lea-Coulson and other models which describe the Luria-Delbrück ﬂuctuation test. A key ﬁnding is that mutant number distributions arising from a generalised Lea-Coulson model for which normal cell growth is non-decreasing are unimodal. An integral criterion is given which separates the cases of a mode at the origin, or not.


Introduction
The Luria-Delbrück fluctuation assay is widely used to estimate mutation rates of micro organisms such as bacterial cells. In very broad outline, several test tubes containing a liquid nutrient medium are seeded with the same number N 0 of normal-type cells. These cells multiply by binary fission attaining the number N t by time t. At this time the contents of the test tube are 'plated' onto a solid substrate which is (almost) immediately lethal for the normal cells. Some cells may have mutated during the growing phase into a resistant type. Under ideal conditions they will form visible colonies on the lethal substrate. Counting these colony numbers provides data which is used to determine the rate of mutation intrinsic to the organism of interest.
Exactly how the data are used for this determination depends on the mathematical model chosen to describe the dynamics of the situation. Various choices are available and we refer to [1,2] for reviews and references. The Lea and Coulson [3] model and its subsequent tweaks is the most widely used of those available. In its simplest form it assumes the following occurs within each test tube.
(i) Normal cell numbers increase exponentially fast: N t = N 0 e νt . (ii) Mutation occurs randomly at a rate proportional to N t . Specifically, there is a mutation rate r (per unit time per bacterial cell) such that a mutation event occurs in the interval (t, t + dt) with probability rN t dt + o(dt), i.e., a normal cell converts to a resistant type.
There is no mutation with probability 1 − rN t dt + o(dt). (iii) Mutation events create mutant clones which grows independently of each other according to a linear birth process with split rate µ, i.e., a binary splitting or Yule process. The relative growth rates of normal to mutant cells is denoted by γ = ν/µ. These assumptions give rise to probability distributions for the total number M t of mutants at time t. The model implies that these distributions are infinitely divisible (abbreviated infdiv), i.e., compound Poisson distributions [4]. Our aim in this paper is to investigate the presence of deeper infdiv properties of mutant number distributions. More specifically, are they generalised negative-binomial convolutions (GNBC's)? The answer is interesting in its own right, but a positive answer gives structural insight, in particular that such a distribution is unimodal and it provides criteria which determine if the associated probability mass function is non-increasing or it has a positive mode.
Definitions and basic properties of positive infdiv distributions are reviewed in §2. Useful subclasses of infdiv distributions are characterised by analytical properties of the density (x) (defined for x > 0) of the Lévy measure (c.f. (1)). The subclass of selfdecomposable (SD) distributions is defined by requiring that x (x) is non-increasing. This class is significant because its members have a unimodal density function. A corresponding discrete version is defined, and they are unimodal too. In addition, precise criteria exist which separate the cases of a non-increasing mass function (i.e., a mode at the origin), or the smallest mode is positive. These notions are applied to models where mutant clones grow as deterministic integer-valued functions.
The Lea-Coulson model described above can be generalised to allow normal cell growth to be an arbitrary positive-valued function of time. The resulting mutant number distribution is a mixture of Poisson distributions where the mixing distribution is a continuous infdiv distribution with the special property that (x) is completely monotone. Distributions having this property comprise the so-called Bondesson (BO) class. This notion is introduced in §3, along with its discrete version. Details are provided for the balanced (γ = 1) generalised Lea-Coulson model in which normal cell lines grow according to the logistic (Pearl-Reed) population model.
The generalised gamma convolution (GGC) class of infdiv distributions comprise the subset of BO distributions for which the product x (x) is completely monotone. This implies the inclusion GGC⊂SD, and hence GGC's are unimodal. Poisson mixtures in which the mixing distribution is a GGC comprise the class of generalised negative-binomial convolutions (GNBC's), and they too are unimodal. Relevant definitions and properties are introduced in §4 where it is shown that a mutant number distribution arising from a generalised Lea-Coulson model in which normal cell growth is non-decreasing is a GNBC. This of course applies to the standard Lea-Coulson model as described above, and details are presented in §4 together with precise criteria concerning the modal behaviour of the mutant number distribution; see Theorem 5(a). The section ends with a discussion of shapes of mutant number distributions selected by different estimation methodologies applied to experimental data and also the preservation of the GNBC property when plating efficiency is an issue.
It is often observed that mutations occur during the time of division of a normal cell. This contingency is addressed by branching process descriptions of the Luria-Delbrück set-up. Some details are provided in §5 for the two most common models, those due to Haldane and Bartlett. Mutant number distributions for the Haldane model are not infdiv, whereas they are infdiv for the Bartlett model. However, rather less can be ascertained about fine infdiv properties for this model.
Finally, in §6 we determine infdiv and modal properties of mutant number distributions arising from alternative models discussed by Kepler and Oprea [5], Angerer [6] and Stewart et al. [7]. Some notation may have different definitions in different sections, but no confusion should arise.

Infdiv Distributions and Deterministic Mutant Growth
In this section, we shall review necessary basic ideas of infinite divisibility and selfdecomposability and explore their (limited) applicability to the Lea-Coulson and Armitage models in which mutant numbers are assumed to increase deterministically.
Let X be a non-negative random variable with distribution function (DF) F(x) and Laplace-Stieltjes transform (LST) If F(x) has a probability density function (pdf) f (x), then F(ζ) = ∞ 0 e −ζx f (x)dx. Denote the left-extremity of F by a = inf{x : F(x) > 0} and observe that a ≥ 0. Thus, F(x) = 0 if x < a and F(x) > 0 if x > a. This quantity can be computed from the LST according to Each of the quantities X, F and F are called infdiv if, for each t > 0, the function is the LST of a probability distribution. This implies that, for any positive integer n, X can be expressed as a sum of random variables, X = ∑ n j=1 X j,n , where the summands are independent and they have the distribution determined by F(ζ) 1/n . This encapsulates the idea of infinite divisibility. It is the case that the sum of independent infdiv random variables is itself infdiv.
An infdiv LST has a special canonical form where c(ζ) is called the Laplace exponent and λ(·) is a measure, called the Lévy measure, which satisfies the conditions This means that λ assigns a zero mass to the origin, it may assign infinite mass to any small interval (0, ) but it integrates x at the origin, and it assigns a finite mass to infinite intervals ( , ∞); here is an arbitrary positive number. Functions having the form (1) are called Bernstein functions-see [8], the standard reference. Differentiation of c(ζ) shows that Many common distributions are infdiv: gamma, Pareto and log-normal, to mention a few. For us the most important is the gamma family. We say that the random variable γ(σ) has the standard gamma distribution with shape parameter σ > 0 if its pdf is g σ ( denotes the gamma function (due to Euler); see [9]. The gamma pdf is decreasing in (0, ∞) if σ ≤ 1 and it has a single positive mode at x = σ − 1 if σ > 1. The corresponding LST is (1 + ζ) −σ , equivalently, c(ζ) = σ log(1 + ζ). We stress that infdiv laws can be multi-modal. Remark 1. In many instances a = 0 but we will need the additional generality for subsequent key definitions.
is a distribution function and the Laplace exponent (1) can be written as with the interpretation that where the J i are independent with DF G and N Λ is independent of the summands and it has the Poisson distribution with (rate) parameter Λ (and denoted by Poisson(Λ)). Thus, X is represented as a (Poisson) random sum of independent jumps J i and it is said to have a compound Poisson distribution. Conversely, any positive infdiv distribution can be realised as the limit of a sequence of compound Poisson distributions.
An important sub-class of infdiv distributions is the class of self-decomposable (SD) distributions. This notion can be given three equivalent definitions but we concern ourselves with the two which fit with our theme. The definition which explains the terminology is that X has a SD distribution if it has the autoregressive representation that, for any constant c ∈ (0, 1), there is a random variable X c independent of X such that This says that if X is scaled down to cX, then the distribution of X can be recovered by adding an independent 'error' X c . Thus, the right-hand side represents the 'selfdecomposition' of X. This definition can be expressed in terms of the LST F(ζ) of X as the assertion that X has a SD distribution if, for each c ∈ (0, 1), the quotient F(ζ)/ F(cζ) is completely monotone, and hence is the LST of a random variable, X c say. It can be proved that a SD distribution is absolutely continuous and infdiv. (In addition, the 'error' term X c is infdiv.) The Lévy measure takes a special form which characterises SD distributions and which sometimes is adopted as the definition of this concept. We shall do likewise with the following formal definition refining (1).

Definition 1.
An infdiv distribution is SD if its Lévy measure λ has a density, where k(x) is non-increasing in (0, ∞). The regularity properties of λ then require that It follows that E(X) = a + ∞ 0 k(x)dx.

Example 1. The integral representation
(just differentiate each side) implies that the gamma(σ) distribution is SD with k(x) = σe −x .
The following fact is important.
(a) Sums of independent SD random variables are SD. (b) If F is the DF of a SD distribution, then it has a pdf f which solves the integral equation This pdf is unimodal, and if k(0+) ≤ 1, then it is non-increasing with a mode at zero. If k(0+) > 1, then f is bounded. In addition, with b := sup{x > 0 : See [10] (pp. 408, 409) for the modality assertions, and more.

Remark 2.
The integral Equation (6) has a wider applicability than is indicated by Fact 1. Specifically, if f (x) is a pdf for which there exists a function k(x) such that (6) holds, then f is infdiv iff k(x) ≥ 0. See [11] (p. 95) for an even more general account.
Since members of the class of SD distributions have an absolutely continuous DF, we may wonder about discrete analogues of this concept. Suppose that X is infdiv and it can take only non-negative integer values, i.e., it is discrete infdiv. Then it necessarily has a compound Poisson distribution with positive integer jumps. J i . Denoting the PGF of the jump distribution by h(s) = E(s J i ), the general form (2) becomes where the notation on the left-hand side anticipates the application of these concepts to mutant number distributions. Here we understand that h(0) = 0, i.e., there are no zero-sized jumps.
Writing the jump PGF as h(s) = ∑ ∞ j=1 h j s j , then setting s = e −ζ in (7) and comparing the result with (1) (with a = 0) makes it clear that the Lévy measure inherent in (7) assigns mass λ j = Λh j to integers j = 1, 2, . . . . Hence the total mass of the Lévy measure is It is often convenient to express the PGF M in the form noting then that logarithmic differentiation of (7)/(8) yields M (s) = M(s)R(s) with r j s j and r j := (j + 1)λ j+1 .
We thus obtain the discrete analogue of (6), Remark 3. The sequence (r j ) is called the canonical sequence, or r-sequence, of the infdiv distribution (p j ). In fact, for any discrete distribution there is a sequence (r j ) such that (10) holds. An essential fact here is a theorem of Katti [12] asserting that (p j ) is infdiv iff its r-sequence is non-negative. See [11] (p. 36). This result has subsequently been 're-discoverd', e.g., [13] and [14] (p. 174).
Many specific discrete distributions discussed in this paper arise as Poisson mixtures where the mixing distribution is infdiv, i.e., where X is infdiv with Laplace exponent (1). Hence Thus the shift term aζ in (1) induces a Poisson(a) component in the discrete mixture. Manipulation of the integral will show that M(s) has the compound Poisson form (7) with A result of Holgate [15] asserts that if the mixing distribution is unimodal (infdiv, or not), then the Poisson mixture is unimodal. The next definition is suggested by Definition 1.

Definition 2.
The discrete compound distribution (p j : j ≥ 0) is called discrete self-decomposable (DSD) if its r-sequence (r j : j ≥ 0) is non-increasing.
The auto-regressive characterisation (4) of (continuous) SD distributions has the following analogue. The characterisation (4) of SD distributions involves multiplying a random variable by the constant c to give a product smaller than X. If X is discrete, then this cannot be done in a way which gives an integer-valued product. Binomial thinning is an analogue which addresses this issue: Define a 'discrete product' as follows. Let p ∈ [0, 1] and where the summands are independent with the Bernoulli(p) distribution and they are independent of X. Thus, E[p X] = pE(X), and the PGF of the product is This product concept is due to the authors of [11]; see p. 495 for the original reference. Definition 3. The discrete random variable X has a DSD distribution if, for each p ∈ (0, 1), there is a discrete random variable X p such that where the summands on the right-hand side are independent. Equivalently, the quantity M(s)/M(1 − p + ps) is a PGF.

Remark 4.
Fact 2 imparts useful qualitative information about the general shape of the mass function of a DSD distribution. If p 0 > p 1 , then p 0 > p j ≥ p j+1 for all j ≥ 1; the mass function is non-increasing. If p 0 < p 1 , then the modal value is positive and it may not be unique. See Discussion 1.
We now consider two models in which normal cells and mutation occur as in §1 and in which mutant clones grow deterministically with sizes having integer values. The first such model was introduced by Lea and Coulson [3] who derived some approximate results for it. Armitage [16] gave it a more careful consideration. More detail is provided by Crump and Hoel [17], who identify it as their D/D 1 model. The survey [1] names it the discretised Luria-Delbrück formulation and the treatment there probably is the most detailed.
Zheng's term captures the central conception that at time t after its formation, the size of a mutant clone is where [·] denotes the 'integer part of'. He shows that the PGF of M t is given by where m = m(t) = (rN 0 /ν) e νt − 1 , θ = θ(t) = (rN 0 /ν)e νt and γ = ν/µ. Theorem 1. The mutant number distribution is DSD, hence unimodal, if (a) K = 1 (equivalently, νt < γ log 2), in which case its mass function is non-increasing iff m ≤ 1; i.e., νt ≤ log(1 + ν/rN 0 ); or if (b) K ≥ 2 and γ ≥ γ * ≈ 0.3663, in which case its mass function is non- Proof. It follows from the definition of K that e µt = K + δ, where δ ∈ [0, 1) is the fractional part of e µt . Hence e −νt = (K + δ) −γ . Substituting into (13) and with reference to (8), a differentiation yields the evaluations If K = 1, then the sum term in (13) vanishes and M t has a Poisson distribution with parameter m and Assertion (a) is known.
If 0 < γ < 1, then this representation of ψ γ is not informative because now the first factor is increasing. Instead, computation of −ψ γ (x) and letting u = x −1 ∈ (0, 1] will show that the sign of −ψ γ coincides with that of Both of σ ± are concave-increasing and hence they can achieve equality in (0, 1] for at most one value of u. Numerical calculation shows that σ + (1) = σ − (1) if γ = γ * specified in the assertion, and that σ Hence the r-sequence is non-increasing if γ ≥ γ * , and Assertion (b) follows from Fact 2.
The case γ ≥ 1 covers the biologically more likely situation in which mutant clones grow no more quickly than normal clones. Theorem 1 fails if γ is sufficiently close to zero. Numerical calculation shows that there is a critical value γ 0 ≈ 0.284 such that r 0 < r 1 (resp. >) if γ < γ 0 (resp. >). In other words, the modal value of the r-sequence jumps from zero to unity at γ = γ 0 . There is a similar jump from 1 to 2 at a critical value γ = γ 1 ≈ 0.179. These outcomes suggest the existence of a sequence of critical values γ i ↓ 0 as i ↑ ∞ at which the modal value of the r-sequence jumps from i to i + 1. In addition, it suggests that The second model we consider derives its deterministic growth character from assuming that mutant cells have a fixed lifetime of duration L at the end of which they divide.
Thus, a clone has size 2 j during the interval [jL, (j + 1)L) since its inception. In order that mutant clones achieve splitting rate µ, we choose L such that Lµ = log 2, i.e., Lν = γ log 2.
This model with γ = 1 was introduced in [17] where it is designated as the D/D 2 model. The expression (11) in this reference for the mutant number PGF is valid for γ > 0 and, with our notation, it is where m and θ are the above time-dependent parameters and now K = [t/L]. We have the following result.
Theorem 2. The mutant number distribution specified by (15) is DSD, and hence unimodal if (a) t < L, in which case the mass function is non-increasing iff m ≤ 1, i.e., Proof. If t < L, then K = 0 and no mutant has reproduced. Thus, M t equals the number of mutations during (0, t] and hence it has a Poisson(m) distribution. Assertion (a) follows.

Bondesson Classes and the Generalised Lea-Coulson Model
In this section, we introduce the first of two special classes of infdiv distributions. The history of these notions is that the Swedish actuary/mathematician Olaf Thorin introduced in 1977/78 distributions now called Generalised Gamma Convolutions (GGC's) with the specific purpose of proving that Pareto and lognormal distributions are infdiv. Subsequently many other distributions conjectured to be infdiv have been proved to be so by showing they are GGC's. A nett benefit of this is that GGC's are SD and hence unimodal. It follows then from Holgate's theorem that Poisson mixtures of GGC's are unimodal too. Lennart Bondesson introduced in 1981 the larger class of infdiv distributions which we review in this section. Detailed accounts of these topics are [18], [11] (Chapter VI) and [8] (Chapters 6-9).
We begin as follows. Let G be a DF on [0, ∞) and define a mixture of exponential distributions by Clearly f is a pdf and the corresponding LST is

Definition 4.
A function F is the DF of a mixture of exponential distributions (written F ∈ ME) if where α ∈ [0, 1] and G is a DF on (0, ∞).
where ε has an exponential distribution and where b(y) is a (measurable) function on (0, ∞) satisfying It follows from Example 1 that the Lévy density of the gamma(σ) distribution is (x) = (σ/x)e −x and, in particular, that it is completely monotone. This motivates the following definition of the class BO of distributions named after Lennart Bondesson.

Definition 5.
An infdiv DF F belongs to the Bondesson class (written F ∈ BO) if its Lévy measure has a completely monotone density, where B is a measure (the Bondesson measure) satisfying where B is a Bondesson measure. (b) The class BO is the smallest set of distributions containing ME and which is closed under convolution and weak limits.
There is a clear similarity of the cumulant functions (16) and (19) with a = 0. This is not mere coincidence.

Fact 5.
If F ∈ BO, then F ∈ ME iff a = 0 and B(dy) = b(y)dy, where b satsfies (17). Definition 6. The discrete random variable X has a geometric mixture distribution if its PGF has the form where Π is a random variable satisfying P(0 < Π < 1) = 1.
is independent of the random variable ε which has a unit exponential distribution, then it follows from the mixture representation of the geometric distribution that The product εC is infdiv, hence any geometric mixture is compound-Poisson.
We now introduce a discrete version of BO; the class BOP of Poisson mixtures with mixing distribution in BO. We will see that mutant number distributions arising from a generalisation of the Lea-Coulson model (below) and from the Bartlett model ( §5) live in BOP.
The following fact arises fairly readily from (19) and Definition 7.
Remark 5. The substitution u = (1 + y) −1 will make clear that (λ j+1 ) really is a Hausdorff moment sequence. For example, if B has a density b(y), then where, in general, δ ρ denotes the measure which assigns unit mass to the real number ρ and zero mass to any interval not containing ρ. The representation asserted in Fact 6 often is more convenient for our purposes.

Remark 6.
In the most general situation, the fact that jump probabilities h j of a compound Poisson distribution comprise a non-increasing sequence implies little about the modal properties of (p j ). For example, if X has the Poisson(h 1 ) and Y/2 the Poisson(h 2 ) distributions, respectively, and X and Y are independent, then X + Y has at least two modes, one at j = 0 and the other at By definition a generalised Lea-Coulson model admits any (measurable) deterministic growth function N(t) of normal type cells. Replacing the exponential form N t with N(t) in the specification of §1 yields a compound Poisson distribution for mutant numbers M t whose Lévy masses are and These outcomes are well-known and they follow from the order statistics property of Poisson processes. See [19] for what seems the earliest and most general formulation. A later independent account specifically for the Luria-Delbrück context is in [17], and the model is reviewed in [1]. This generalised Lea-Coulson model can also be regarded as a branching process with inhomogeneous immigration. The branching component comprises the independently growing mutant clone birth processes and immigrants comprise the inhomogeneous Poisson process of mutations. See [20] for a review of this topic. We have the following general result.
Theorem 3. Let t > 0. The mutant number distribution of the generalised Lea-Coulson model is a BOP distribution whose Bondesson measure has the density b(y, where φ = 1 − e −µt and H(x) = x + is the Heaviside unit step function.
Proof. Just make the substitution (21) and refer to Fact 6 to obtain the desired moment representation, The resulting infinite integral does converge because it equals the integral (21). Alternatively, observe that b(∞−, t) = (r/µ)N(t) and lim y→φ −1 −1 b(y, t) = (r/µ)N(0), implying that the regularity conditions in Definition 5 always are satisfied.
For computational purposes it is more convenient to shift the integration variable in Fact 6 to obtain and the corresponding Lévy density where Remark 7. Substituting, again, u = y −1 in (22) and (24) gives the 'explicit' moment representation Hence the representing measure for any mutation number distribution derived from a generalised Lea-Coulson model has the time-dependent support [0 This moment relation yields the fundamental relations and hence and Example 2. Suppose that normal cells increase in number as a logistic growth model with carrying capacity K > 0. Thus, whose well known solution is . (29) Hence, for the balanced case, µ = ν, some manipulation yields We assume that N 0 < K, implying that N(t) < K, B < 1 and lim K→∞ BK = N 0 e νt .
Substitution into (27) leads to the explicit form where E 1 (x) = ∞ x y −1 e −y dy is the exponential integral; see [9] (# 6.2.1). Define θ K = (r/µ)BK. The integrand of (26) resolves into partial fractions: It follows that and We obtain expressions for the Poisson rates as follows. Writing where, as usual, ∑ 0 i=1 (·) = 0. It follows that The power series expansion of the logarithm term yields the form Letting K → ∞, recalling that B → 0 and noting that θ K → (r/ν)N 0 e νt recovers the balanced Lea-Coulson model which we will consider in more detail in the next section.

Thorin Classes and the Lea-Coulson Model
We now introduce the above mentioned GGC class of infdiv distributions which are pertinent to a significant subclass of generalised Lea-Coulson models. We motivate the general definition by observing that, given independent gamma random variables γ(σ i ) (i = 1, . . . , n) and constants c i > 0, it follows from (5) that the Laplace exponent of the sum where U n is a measure which assigns mass σ i to the point c i . It follows from Fact 1(a) that X is SD. The SD class is closed under limits in distribution so, taking the informal limit n → ∞ in (32) yields a putative limiting Laplace exponent This does specify a SD distribution for any measure U on (0, ∞) satisfying Definition 8. A distribution whose Laplace exponent has the form (33) where a ≥ 0 and U is a measure on (0, ∞) subject to (34) is called a generalised gamma convolution (GGC). A function of the form (33) is called a Thorin Bernstein function. An equivalent specification is that the class of GGC's is the smallest which contains scaled gamma distributions and is closed under convolution and weak limits.
The representing measure U in (33) is often called the Thorin measure and we define the Thorin distribution function T(y) = U((0, y]). We motivate a discrete version of GGC s by observing that the best known case of a Poisson mixture (12) is where X = γ(σ)/c, c a positive scaling constant, giving where p = (1 + c) −1 ∈ (0, 1). Hence this gamma-mixed Poisson distribution is the negative binomial distribution with parameters p and σ, denoted NB(p, σ). The case σ = 1 of course is a geometric distribution whose mixing distribution is an exponential one. The following definition extends this idea.

Definition 9.
A Poisson mixture distribution is a generalised negative-binomial convolution (GNBC) if the distribution of the mixing random variable X is a GGC as defined above.
A calculation using Fact 7 gives where V is a right-continuous function on (0, 1) such that V(1−) = 0, (b) The r-sequence (c.f. (9)) is a Hausdorff moment sequence, Conversely, if the r-sequence of a DID distribution has this moment representation, then it is a GNBC. (c) The GNBC class is the smallest class of discrete distributions which contains negative-binomial distributions and is closed under convolution and weak limits. (d) A GNBC is discrete unimodal and its mass function is non-decreasing iff λ 1 = a + 1 0 udV(u) ≤ 1.

Remark 8.
Since the shift constant a in (33) induces a Poisson(a) component in (35), the leftextremity of a GNBC always is zero. Assertion (d) follows from Fact 7(b) and Holgate's theorem [15], and then Fact 2 observing that r 0 = λ 1 .
The following fact gives a canonical representation for a mixture of geometric distributions and a condition that it be a GNBC; [11] (pp. 381, 390).
Referring to (35), we will later need a general relation between the function V and the Thorin measure U of the mixing GGC distribution. The following result achieves this in terms of the Thorin distribution function T(x).

Proof. The integral in (33) can be written as the Stieltjes integral
It follows from the first member of (34) that for any x, ∈ (0, 1) we can choose δ ∈ (0, 1) such that log y −1 dT(y) < .
Next, it follows from the second member of (34) that there exists x > 1 such that if x > x , then implying that lim x→∞ x −1 T(x) = 0.
Observing that the integrand in (37) is asymptotically proportional to log x −1 as x → 0, and to x −1 as x → ∞, it follows from an integration by parts that In a similar manner, it follows from (35) with a = 0 that the PGF of the corresponding GNBC is The left-hand side equals c(1 − s) and a computation shows that − log M(1 − ζ) reduces to a Stieltjes integral as above with T as asserted.
Recall the expression (22) for the Lévy masses λ j pertaining to the generalised Lea-Coulson model. A very natural condition on the growth function N(t) of normal cells implies that mutation number distributions are GNBC's.

Theorem 5. Assume that the normal cell growth function N(t) is non-decreasing. Fix t > 0. Then: (a) The distribution of M t is a GNBC and hence unimodal. Its mass function is non-increasing iff
(d) The mutant number distribution is a geometric mixture iff rN(t) ≤ µ.

Remark 9.
It follows from (20) and the hypothesis of Theorem 5 that λ 1 is an increasing function of t. Clearly λ ≤ 1 if t is sufficiently small in which case the mutant number mass function will be non-increasing. It attains a positive maximum value if λ 1 eventually exceeds unity.
The logistic differential Equation (28) implies that if N 0 < K, then its solution is strictly increasing. It follows that the corresponding mutant number distribution is a GNBC. However, except for the balanced case γ = 1 it does not seem that the integrals (26) and (27) can be evaluated in any insightful way. In the balanced case we now know that the Lévy density (30) is such that x (x) is completely monotone. The following direct demonstration of this fact yields its Thorin function T(y, t).
Integration by parts shows that The substitution y = (φ −1 − B)v − (1 − B) exhibits x (x) as the sum of two completely monotone functions: Comparing this with (38) we see that Thus the Thorin measure has a discrete component -a point mass at y = φ −1 − 1 and its support is independent of K.
In the remainder of this section we restrict consideration to the Lea-Coulson [3] model described in §1 and give a self-contained treatment starting from (26). Taking N(t) = N 0 e νt we thus obtain where In the sequel we usually suppress the time dependence, thus regarding the distributions determined by (40) as a parametric family determined by (θ, φ, γ) where φ ∈ (0, 1) and θ, γ > 0. Expressions equivalent to (40) appear first in [21]. Sometimes [22] is coupled with this reference because, independently, a system of differential equations for the mass function of M t is derived, generalising the system in [3] for the case γ = 1, and deducing a numerical solution scheme. The integral in (40) has no simple evaluation except perhaps for γ = 1, 2, . . . .
In fact, if γ = 1, then evaluation gives the familiar outcome This PGF appears for the first time in [16] (p. 10) as a result of solving the linear firstorder partial differential equation derived in [3]. Zheng [1] denotes the corresponding distribution by LD(θ, φ) where the LD letter designation is chosen to honour the pioneering contribution of Salvador Luria and Max Delbrück. Frequently in laboratory situations the product µt is so large that φ ≈ 1 and it is argued that the form (42) is approximated by This is a PGF as can be deduced from the explicit time-dependent LD(θ(t), φ(t)) distributions by allowing µt → ∞ (implying φ(t) → 1) and r → 0 such that θ(t) → θ ∈ (0, ∞); a kind of Poisson approximation. Zheng [1] (and others before him) name the distribution corresponding to (43) after Lea and Coulson because they derive (43) by using a clever manipulation to solve their partial differential equation. It is denoted by LC(θ) and thus coincides with LD(θ, 1). The solution (42) satisfies M(s, 0) = 1, reflecting the assumption (and laboratory situation) that M 0 = 0. The LC solution does not satisfy this initial condition, but it has an interesting form-invariant character which bears the interpretation that mutant numbers evolve as a non-homogeneous Poisson process. In view of this historical progression, we will designate the full family of distributions corresponding to (40) by LDM(θ, φ, γ).
It is well known that the LC(θ) distribution is qualitatively very different to LD(θ, φ) distributions when φ < 1. The moments of the former are infinite, reflecting the very slow decrease of its right-hand tail. If φ < 1, then all moments are finite and the right-hand tail decays exponentially fast [4]. The following result shows that each LDM distributions is a GNBC and that the just-mentioned differences are reflected in the representing measures of the mixing GGC. Here, and below, recall that H(x) denotes the Heaviside unit-step function, i.e., the DF of the degenerate distribution allocating unit mass to the origin. Just below, and later, we will encounter the second confluent hyper-geometric function, where a > 0 and b is real. Observe that this function is completely monotone; [9] (Chapter 13). Theorem 6. If γ, θ > 0 and φ ∈ (0, 1], then the LDM(θ, φ, γ) distribution has the following properties.
(a) It is a GNBC, hence unimodal. Its mass function is non-increasing iff In particular, the LDM(θ, φ, γ) distribution is a geometric mixture iff θ ≤ 1.
(c) The GGC mixing distribution has the Thorin distribution function The Lévy measure of the GGC mixing distribution has a density which has the following explicit forms: Remark 10. The above-mentioned difference between the cases φ = 1 and φ < 1 are manifested in the fact that the representing functions V and T are continuous with supports coinciding with their domains iff φ = 1. Indeed, if φ < 1, then −V(u) decreases from θ at u = 0 to θ(1 − φ) γ at φ− and it jumps to zero at u = φ. Note that (48) results by letting K → ∞ in (30).

Proof. (a) Comparing (6) with (40), a differentiation gives
Hence This exhibits the desired Hausdorff moment form with the measure This implies the first assertion, and the second follows by evaluating r 0 = R(0) and appealing to Fact 2. Observe that the measure V has a discrete component which vanishes when φ = 1. (b) Integrating (49) and simplifying the result leads to where C is the constant of integration. The condition V(1−) = 0 implies that C = −θ, whence (45). (c) The evaluation (46) comes directly from Theorem 1 and (45).
(d) Recall that the Lévy density (x) = x −1 k(x) exists and, with no parameter restriction, The right-hand side integral is an 'incomplete' confluent hypergeometric type of integral.

Remark 11.
Reverting to the time-dependent form of parameters, it follows from Theorem 4 that being a geometric mixture and the nature of modality are time-dependent properties, whereas, e.g., the SD property of the LDM distributions is a time-independent property. See [10] for this dichotomy.

Discussion 1.
The LDM(θ, φ, γ) family of distributions is most commonly used to fit empirical mutant number distributions. It follows from the criterion (44) that as θ increases from small to large values, the mutant number mass function transitions from decreasing to having a positive mode. If equality folds in (44), then zero and unity are modes. It usually is the case that the estimate of φ is so close to unity that it is chosen to equal unity. In this case the criterion (44) simplifies to In the case of equal fitness of normal and mutant cells, γ = 1 (the LC(θ) model), then the transition from a zero to positive mode can be seen in the first three columns of Table 2 in [3] where, if θ = 2 (denoted by m in this reference), then p 0 = p 1 = 0.1353 > p 2 = 0.1128. The LC(θ) model is fitted to three sets of laboratory data in [22] where θ is estimated as 0.3783, 3.84 and 3.03, respectively. Figures 3-5 in [22] graph the mass functions corresponding to these values.
Finally, to see that real estimated mutant number distributions can exhibit a zero or a positive mode we recall estimates determined in [23] from several experimental data sets for the LDM(θ, 1, γ) distribution. A main objective in [23] is to introduce parameter estimation based on the empirical PGF and compare its performance with maximum likelihood estimation (MLE). Table 1 in [23] presents 95% confidence intervals for θ (denoted there by α) and γ (denoted there by ρ). Assuming that point estimates are the mid-point values of the confidence intervals, Table 1 here exhibits these estimates and it indicates the shapes of the estimated mass functions. There now are several methods of estimating mutation model parameters and a question of interest is that if several methods are applied to a given set of data, will they be consistent as to the shape of the mutation number distributions they select? Published studies indicate that different methods can give quite different estimates, but they usually are, but not always, consistent in regard to the selected distribution shape. We mention two comparative studies for the LC(θ) model.
Five estimation methods are compared using four data sets in [28] (where m is used for θ). Table 2 in [28] shows broad consistency in shape selection for Experiments A-C, with the first two indicating a zero mode and the third a positive mode. The p 0 -method was not applied/applicable to the Experiment D data. Two of the four estimated values resulted in a mode at zero, and the other two a positive mode. Table 5 in [26] compares seven estimation methods using seven sets of experimental data. Estimates of θ (m in [26]) are quite variable across estimation methods, but selected shapes are broadly consistent. In fact, after adjusting the Luria-Delbrück mean method by eliminating large jackpots, all methods were consistent in five of the seven data sets. In the two other cases all methods except the Drake median method gave estimated θ < 2, and the Drake estimate was a little over two; 2.07 for Experiment 2 and 2.08 for Experiment 6. In these cases the modal value is unity; p 0 = 0.126, p 1 = 0.131 and p 2 = 0.111 if θ = 2.07, and p 0 = 0.125, p 1 = 0.130 and p 2 = 0.111 if θ = 2.08. Finally, a zero mode was found for five of the seven data sets.
These investigations do provide confidence that, although different methodologies can show rather different parameter estimates, they in fact are broadly consistent with respect to shape selection.
We end this section with some remarks about plating efficiency. This term refers to the possibility that, upon plating, a mutant cell fails to establish a colony. This aspect frequently is modelled by assuming that plated mutants independently establish colonies each with a probability p ∈ (0, 1]. In other words, successful establishment is modelled by binomial thinning -if M(s) is the PGF of the number of plated mutants, then the PGF of the number of established colonies is M e (s) = M(1 − p + ps). A very convenient result asserts that binomial thinning preserves the GNBC property. Specifically, if q = 1 − p, then we obtain from (35) that where dV p (w) = dV(u) and u = w p + qw .
In particular, if these measures have densities v p (·) and v(·), respectively, then v p (w) = p (p + qw) 2 v w p + qw .

Branching Process Models
The normal population is depleted by one cell each time a mutation occurs. The Lea-Coulson model does not directly account for this. One argument is that in real situations N t M t so, this contingency can be neglected. Another response is to replace the parameter ν with ν − r, thus adjusting for a diminished average normal population growth rate.
Branching process models do take direct account of the normal population diminution due to mutation. A discrete-time model was propounded (no later than 1946) by J.B.S. Haldane. See Zheng [29] for an account and references. Haldane's model counts population sizes generation by generation. Cell numbers increase by binary division and hence the total size (normals plus mutants) of the nth generation cannot exceed N 0 2 n . Consequently the distribution of M n , the size of the nth mutant generation, cannot be infdiv. There is a Poisson type of limit theorem [30] resulting in a limiting compound Poisson distribution (and hence infdiv) and whose jump distribution has the PGF h(s) = ∑ j≥0 2 −j−1 s 2 j , a gap series, and hence this limiting distribution is not DSD.
Instead, we shall consider the continuous-time version of Haldane's model. This model is a two-type linear birth process apparently formulated by M.S. Bartlett around 1951/2. It is mentioned for the first time in [16] (p. 37) with details appearing in the first edition of [31] (p. 132) published in 1955. See [32] for a detailed account and earlier references.
The balanced version of the model assumes that normal and mutant types divide into two cells during the interval (t, t + dt) with probability µdt + o(dt), independently of previous history. Mutants breed true, but a dividing normal cell has probability p of producing one mutant and one normal cell and probability 1 − p of producing two normal cells.
The PGF of M t is where φ = 1 − e −µt (as above), but again we suppress the dependence on time t in our notation. Zheng [32] (with more detail in [33]) notes a Poisson type of limit in which φ → 1 (i.e., t → ∞) and p → 0 such that resulting in the limiting PGF (51) The following result gathers infdiv properties of the Bartlett distributions, however, it is deficient in NOT concluding that they are GNBC's. Referring to (50), the term in square brackets can be written as [1 + α(1 − h(s))] −1 , where α = pφ/(1 − φ) and say. We show below that h(s) is a PGF. It follows that in (50), the integer N 0 can be replaced with a positive-valued parameter, σ say. Thus i.e., M(s) is the PGF of a gamma mixture of discrete infdiv distributions. We shall denote members of the resulting Bartlett family of distributions by B(φ, p, σ). The next result shows that a Bartlett distribution is a gamma mixture of GNBC's.

Theorem 7.
Let Λ > 0 and h(s) be as defined in (53). The distribution whose PGF is P(s) = exp[−Λ(1 − h(s))] is a GNBC whose mixing GGC has the Lévy density The corresponding Thorin distribution function is Proof. We begin by showing h(s) is a PGF. Writing c(s) = ∑ j≥0 c j s j and referring to (53) we see that c 0 = 1 and where B(·, ·) is the beta function. We thus have the explicit representation , by virtue of a reflection formula for gamma functions. It thus follows from the usual integral representation for beta functions that Hence h is a PGF, as asserted above.
Making the substitution u = (1 + y) −1 and comparing the outcome with (20) we see that the Poisson intensity sequence (λ j+1 : j ≥ 0) is a Hausdorff moment sequence and that the Bondesson measure has support [φ −1 − 1, ∞) and density b(y) = (Λ/K)φ p−1 (1 + y) p · y 1 + y 1 − 1 This and Fact 7(c) imply (55). Referring to (18), we have (x) = (Λφ p−1 /K) η(x) where The second equality above follows from the substitution y = z + φ −1 − 1 and the final form follows from evaluating the subtracted integral term in the penultimate line using the substitution z = y/φ to obtain We thus have obtain a final outcome and it follows from its construction that is completely monotone. Hence P(s) is the PGF of a BOP distribution. Furthermore, this exhibits (x) as the difference of completely monotone functions and we need to find a different representation to be able to conclude that x (x) is completely monotone. 1 − p, ξ). Integration by parts of the right-hand side integral leads to Substitution into (56) yields (54), as asserted. It is clear now that x (x) is completely monotone, and hence that P(s) is the PGF of a GNBC. It follows from Theorem 7 that this involves the composition of two Thorin Bernstein functions. However, the class of such functions is not closed under composition and hence we cannot conclude that a Bartlett distribution is a GNBC. On the other hand, the components of this composition are complete Bernstein functions and this class is closed under composition. See [8] (pp. 112 and 94), respectively. Hence we can conclude that Bartlett distributions of mutant numbers belong to BOP.
Similarly, the Zheng PGF (51) is that of a gamma mixture of Lea-Coulson distributions. Hence a corresponding analogue of Theorem 7 in essence is Theorem 6(d).

Some Other Mutant Number Distributions
The total population size n t = N t + M t for the above (balanced) Bartlett model comprises a linear birth process with splitting rate µ. Thus, the embedded jump chain is the deterministic process which jumps by unity at each cell division. Angerer [6] and Kepler and Oprea [5] independently and almost simultaneously proposed a discrete-time model for mutant numbers M n immediately following successive divisions at which n t takes values n = N 0 , N 0 + 1, · · · . Thus, M 0 = k if n = N 0 + k, and clearly M n ≤ n − N 0 . Their precise specifications differ in some details but, as in §5, a dividing normal cell produces one normal and one mutant with probability p. Angerer mentions back mutation but does not pursue that issue, instead he allows for mutation rates to depend on n and he provides a very careful and exact treatment of their models. Kepler and Oprea include the possibility of back mutation. Taking account of these differences, their fundamental difference equations relating the distributions of M n and M n+1 , Equation (1) in both references, are the same.
In a more detail, Kepler and Oprea [5] assume a dividing mutant produces two mutants with probability 1 − q and one cell of each type with probability q 1. With no detail given, after they 'pass to a continuum representation form', they assert that the PGF M(s, n) of M n is given by − log M(s, n) = pn(1 − s) Let δ = 1 − p − q ∈ (−1.1), although the biological context implies that 0 < δ 1. Note that taking δ = 0 yields the Poisson distribution with parameter p(n − N 0 ).
So, assuming that δ ∈ (0, 1), the substitution u = 1 − v δ and then comparing the outcome with (40) shows that M n has the LDM(θ, φ, a) distribution with θ = pn/δ, φ = 1 − (N 0 /n) δ and a = δ −1 > 1, and hence Theorem 6 above applies. Angerer [6] proves several limit theorems for M n as n → ∞ and other constraints hold. For example, the limiting PGF displayed as (29) in [6] shows that the limiting distribution is that of a sum N B + M, where N B has a negative binomial distribution and M a LD distribution and they are independent. Hence the sum is a GNBC. Similarly the limit (32) in [6] is the PGF of a similar sum with N B replaced by a Poisson distributed random variable, again a GNBC.

Proof. (a) We have
Expanding the logarithm term and collecting coefficients of s j , we find that Λ = ϕ(1 − p) and Thus the sequence (λ j+1 : j = 0, 1, . . . ) is a Hausdorff moment sequence, implying membership of BOP.
Recalling that r j = (j + 1)λ j+1 , we have Hence we obtain a moment representation r j = ϑ This function increases in (p, 1) but it has a negative jump at u = p. Hence it is not monotone, implying the second assertion of (a). (b) The second equality of (59) can be expanded as Hence r 0 ≥ r 1 iff p ≤ 1/4, a necessary condition for the SD property. In addition, r j−1 ≤ r j iff The left-hand side of (60) is bounded below by Hence (60) certainly holds if p ≤ j/2(j + 1). The right-hand side is increasing in j and the case j = 1 requires that p ≤ 1/4. So, this condition is sufficient for the SD property.
The final model we shall examine is based on the discretised Luria-Delbrück model as reformulated in [7]. There are three model assumptions: 1.
The probability of a mutation during (t, t + dt) is φ(t)dt, where φ(t) = rN t , but otherwise is arbitrary; 2.
A mutation occurring at time t induces a growing clone of size C t at the time of plating/observation. Define p(j, t) = P(C t = j); and 3.
Mutations are classifies as type j if C t = j. The number of type j mutations in a single culture is denoted by M j , a random variable having a Poisson(λ j ) distribution where and "T is the time after which no observable mutations will occur". Presumably, this could be the time of plating.
In relation to the second assumption, there is an enigmatic assertion that p(j, t) "depends on when the mutation occurs". However, this is the absolute time t according to their direct specification. So, perhaps what is meant that t here means the current lifetime of the clone. We shall adopt this interpretation because it seems best aligned with the third assumption. Thus, M j is the number of type j mutations existing at time T.
Consequently, the number of mutants at time T is M = ∑ j≥1 jM j and, assuming that the M j are independent, which is unstated but implicit in [7], the PGF of the mutant number distribution is where, as above, Λ(s) = ∑ j≥1 λ j s j and Λ = Λ(1). Thus, the computation of M(s) reduces to a determination of φ(t) and p(j, t).
A Luria-Delbrück model with a time and state-dependent mutation rate is specified in [7] (p. 181). Normal cell numbers grow according to N t = N 0 e νt and mutant numbers as where, in the integral, we regard t as a function of n = N t .
Thus, the problem reduces to deciding the form of φ(t). A standing assumption is that φ(t) = rn and dn = νndt and, more specifically, that φ(t)dt = rdn + αndt, where dn and dt are related by and α, P and Q are positive constants. Here, rdn represents a constant mutation rate per cell per generation and αndt a rate per cell per time. These specifications yield and hence Observe that the integral for λ 1 diverges for j = 1. This is handled by computing the rate λ 1 required to achieve a specified value of Λ, although this tactic does represent a deviation from the model formulation in [7]. The above log-term equals log(j + 1) + log(j − 1) − 2 log j, and hence partial summation yields Note that an approximation has been adopted in [7] whereby the zero-valued λ j are replaced by the algebraic values obtained from the integration. It follows that a necessary condition for DSD is r 0 = λ 1 ≥ r 1 = 2λ 2 , i.e., Λ ≥ AN T /6 + B log(32/9) = AN T /6 + 1.2685B.
Theorem 9. The mutant number distribution for the above specification is DSD iff (62) holds, in which case the mass function is non-increasing iff λ 1 ≤ 1.
Proof. If j ≥ 1, then The coefficient of B equals For any v ∈ (0, 1), the integrand decreases as j increases from (1 − 1 2 v) −1 to (1 − v) −1 and the length of the interval of integration decreses too. Hence r j > r j+1 if j ≥ 2.
We know that the sequence of Poisson rates whose terms equal (j(j + 1)) −1 correspond to a GNBC. So a question is whether the sequence of rates − log(1 − j −2 ) (j ≥ 2) together with an admissible value for λ 1 similarly can be associated? We shall see below rhat the answer is No! Referring to (61), if P, Q → ∞ such that Q/P → 1, then the result is the differential equation for logistic growth. Hence (61) itself represents a generalised form of logistic growth. More generally, (61) is a particular case of the relation dn = nL(n)dt, where L(n) is decreasing in n. We choose the following specific form.
Let ε ∈ (0, 1) and Thus ε = 0 gives logistic growth, and if 0 < ε 1, then L(n) has a quadratic profile approximating the linear logistic profile. We compute Evaluation of the integral follows from the substitution v = n/N T and resolving the integrand into partial fraction form. Note that the cases ε = 0 and j ≥ 2 yield the sequence in [7] and that our restriction ε < 1 is required by the context because L(n) is increasing if ε ≥ 1.
Proof. (a) Begin with the following easily checked identity The integrand term in brackets equals 1 0 u c−1 u − √ y − u √ y du.
Thus we obtain a double integral and the integral with respect to y is It follows from c − 1 > 0 that lim u→0 u c−1 / log u = 0 and in addition, lim u→1 (1 − u) 2 / log u = 0. The first equality in (65) follows, and a log-differentiation yields the second equality.
it follows that the integral expression for the Lévy density of the mixing GGC is