Electronic Journal of Differential Equations, Vol. 2021 (2021), No. 56, pp. 1–22. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu AGE-DEPENDENT BRANCHING PROCESSES AND APPLICATIONS TO THE LURIA-DELBRÜCK EXPERIMENT STEPHEN J. MONTGOMERY-SMITH, HESAM OVEYS Abstract. Microbial populations adapt to their environment by acquiring advantageous mutations, but in the early twentieth century, questions about how these organisms acquire mutations arose. The experiment of Salvador Luria and Max Delbrück that won them a Nobel Prize in 1969 confirmed that mutations don’t occur out of necessity, but instead can occur many generations before there is a selective advantage, and thus organisms follow Darwinian evolution instead of Lamarckian. Since then, new areas of research involving microbial evolution have spawned as a result of their experiment. Determining the mutation rate of a cell is one such area. Probability distributions that determine the number of mutants in a large population have been derived by Lea, Coulson, and Haldane. However, not much work has been done when time of cell division is dependent on the cell age, and even less so when cell division is asymmetric, which is the case in most microbial populations. Using probability generating function methods, we rigorously construct a probability distribution for the cell population size given a life-span distribution for both mother and daughter cells, and then determine its asymptotic growth rate. We use this to construct a probability distribution for the number of mutants in a large cell population, which can be used with likelihood methods to estimate the cell mutation rate. 1. Introduction 1.1. Luria-Delbrück experiment. In the early twentieth century, questions about how microorganisms acquire advantageous mutations arose. In 1943, biologists Sal- vador Luria and Max Delbrück conducted an experiment in order to determine whether mutations occurred out of necessity or could occur many generations be- fore there was a selective advantage. This experiment, dubbed the “Luria-Delbrück Experiment,” helped them win a Nobel Prize in 1969 (see [6] and [12]). In their experiment, Luria and Delbrück grew bacteria in a non-selective medium in multiple tubes for a period of time until they all reached a certain cell density. Then, they plated the cells from each tube on different plates of a selective-medium containing a bacterial virus. Cells that showed resistance to the virus had acquired a virus-resistant mutation. If cells evolved according to a post-exposure hypotheses such as Lamarckian evolution, where cells acquire mutations in response to their environment, then the number of mutants in each plate would follow a Poisson 2010 Mathematics Subject Classification. 92D15. Key words and phrases. Probability generating function; fluctuation analysis; asymmetric cell division; Laplace transform. c©2021. This work is licensed under a CC BY 4.0 license. Submitted May 26, 2021. Published June 23, 2021. 1 2 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 process where the mean is equal to the variance, making plates with a large number of mutants highly unlikely. But in the experiment of Luria and Delbrück, there were “jackpots,” meaning there were plates with an unusually large number of surviving cells. The only conclusion they could make was that mutation occurred before the cells were plated on the selective medium containing the virus, and thus Charles Darwin’s theory of natural selection applied to microorganisms. Since then, new areas of research involving microbial evolution have spawned from the Luria-Delbrück Experiment which are still being studied today. Since cells can acquire mutations before there is a selective advantage in their environ- ment, questions about their mutation rate have risen. However, traditional methods involving significance tests can not be used due to the high variability of the data. Instead, we can use likelihood methods, but in order to do so, we need a probabil- ity distribution for the mutant cell population size as a function of the unknown mutation rate. Though unpublished originally, a probability distribution for the number of mu- tants was presented by Haldane (see [11]). However, Haldane’s model had two major issues: it assumed all cells divide synchronously, and the distribution was computationally inefficient. In 1949, Lea and Coulson constructed a generating function for the number of mutants in a large population that had a closed-form solution, so computing the probability coefficients was much more efficient. How- ever, they assumed that all cells grow symmetrically with a life-span distribution that was exponential, which is a very broad assumption for how cells grow (see [5]). In this paper, we look to extend on the ideas of Haldane, Lea, and Coulson and develop a probability distribution for the mutant cell population size where we have control of a cell’s life-span distribution. In addition, we develop a distribution for asymmetric cell division, where a cell divides into a mother and daughter cell with different life-span distributions. 1.2. The problem. Suppose you start with a single cell, and this cell undergoes binary division and divides into two cells. As time passes, the total cell population will grow, but depending on when each cell divides, the population will vary. If a cell divides into two cells identical to itself (symmetrical cell division), then the two children cells will divide similarly to its parent. On the other hand, how will the cell population grow if a cell divides into two cells, where one is identical to itself, but the other is not (asymmetrical cell division)? In addition, suppose during any cell division, a mutation can occur. If this mutation is passed through all the children spawned from the mutated cell, what can we say about the distribution of the mutant cell population size in a sufficiently large cell population? The solution to this problem is given in Theorems 5.8 and 5.9, in the form of probability generating functions. In Section 6 we show how to explicitly compute these probabilities. 2. Probability generating functions A very brief overview of probability generating functions is given to ensure the reader is familiar with the basic properties. The first thing to note is that a prob- ability generating function can only be constructed for discrete random variables outputting non-negative integers. EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 3 Definition 2.1. Let X be a discrete random variable outputting non-negative integer values. The probability generating function (or simply generating function) of X is the function GX(z) := E [ zX ] = ∞∑ k=0 Pr(X = k)zk. (2.1) Probability generating functions are power series with non-negative coefficients such that their sum is 1, so their radius of convergence is always at least 1. Proposition 2.2. If X is a discrete random variable outputting non-negative in- teger values with generating function GX(z), then E[X] = G′X(1). (2.2) Below, we define joint probability generating functions, since we will use them in Section 3. Definition 2.3. Let X and Y be jointly distributed discrete random variables out- putting non-negative integer values. Then the joint probability generating function is GX,Y (x, y) := E[xXyY ] = ∞∑ k=0 ∞∑ j=0 Pr(X = k, Y = j)xkyj . (2.3) Proposition 2.4. If X and Y are independent random variables, both outputting non-negative integer values, then the generating function of X + Y is GX+Y (z) = GX(z)GY (z). (2.4) The following two propositions will be used extensively in Section 5 when con- structing the generating function for the mutant cell population. Proposition 2.5. If X = Xk are independent and identically distributed random variables outputting non-negative integer values, and Z = N∑ k=0 Xk, (2.5) where N is an independent random variable outputting non-negative integer values, then the generating function of Z is GZ(z) = GN (GX(z)). (2.6) Proposition 2.6. If X is a Poisson random variable with parameter λ, then the generating function for X is GX(z) = eλ(z−1). (2.7) For a more in-depth look at generating functions, one can refer to almost any book on probability theory, though we only use the above definitions and properties in this paper. 4 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 3. Cell growth under asymmetric division It is of great interest to biologists to model the growth of a cell under asymmetric division; that is, when it divides into two types: a mother cell and a daughter cell. The mother cell is a copy of the parent cell, while the daughter cell typically takes time to grow into a mother cell. Therefore, the life-span distribution for mother cells is different, and often their life-span is shorter than daughter cells. In this section, we will construct time-dependent generating functions (3.3) for the mother and daughter cell population sizes when we start with exactly one mother cell and when we start with exactly one daughter cell. These generating functions will be in the solution of an integral equation, so we will proceed to show existence and uniqueness of solutions in Section 3.3. We conclude the section with some examples of generating functions using different life-span distributions. In [9], the symmetric case is expounded in full detail before the asymmetric case is described. In this paper, in the interests of brevity, we only explain the asymmetric case, since the symmetric case follows as a special case. 3.1. Preliminaries and assumptions. We will start by stating our assumptions about asymmetric cell division. Assumption 3.1. Cells have the following properties: (1) there are exactly two types of cells: mother cells and daughter cells; (2) all cells are independent of each other, mother cells are identical to other mother cells, and daughter cells are identical to other daughter cells; (3) cell life-span for mother cells and daughter cells are strictly positive, real- valued random variables Tx and Ty, respectively, with distributions P (t) = Pr(Tx ≤ t) Q(t) = Pr(Ty ≤ t) (3.1) respectively; (4) at the end of a cell’s life, both mother cells and daughter cells will divide into one mother cell and one daughter cell. Let Xt and Yt be random variables representing the mother and daughter cell populations at time t ≥ 0, respectively. When we start with exactly one mother cell and no daughter cells, which is the main case of interest, we will have X0 = 1 and Y0 = 0. On the other hand, when we start with exactly one daughter cell and no mother cells, we will have X0 = 0 and Y0 = 1. Note that Xt and Yt are not independent processes. 3.2. Constructing the generating function. To derive the generating function for the cell population at any time t ≥ 0, we will first construct a model in discrete- time, then divide our time increments infinitesimally small, and finally take limits to derive a continuous-time model. Suppose first that time is discrete and the life-span random variables for mother and daughter cells, Tx and Ty, only output positive integer values. For clarity, when we say a cell divides at time t, we mean the cell population has increased by one at time t+ 1. Define f(t, x, y) := E[xXtyYt | X0 = 1, Y0 = 0] g(t, x, y) := E[xXtyYt | X0 = 0, Y0 = 1] (3.2) EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 5 to be the joint generating functions of Xt and Yt, with different initial values. Then it can easily be shown that f(t, x, y) = xPr(Tx > t) + t∑ k=1 f(t− (k + 1), x, y)g(t− (k + 1), x, y) Pr(Tx = k) g(t, x, y) = yPr(Ty > t) + t∑ j=1 f(t− (j + 1), x, y)g(t− (j + 1), x, y) Pr(Ty = j). Now, to get the formulas in continuous time, we take a limit of the discrete time formula as the discrete-time increments converge to zero, and we obtain f(t, x, y) = x(1− P (t)) + ∫ t 0 f(t− τ, x, y)g(t− τ, x, y) dP (τ) g(t, x, y) = y(1−Q(t)) + ∫ t 0 f(t− τ, x, y)g(t− τ, x, y) dQ(τ). (3.3) Since the integrals are Lebesgue-Steiltjes integrals, it is important to emphasize that ∫ t 0 denotes integration over the closed interval [0, t]. 3.3. Existence and uniqueness. A variant of arguments given in [4], Chapter VI, §9 will also work, but we present a slightly different approach. Showing that there exists two unique generating functions f and g that satisfy (3.3) given any life-span distributions P for mother cells and Q for daughter cells will follow by constructing suitable functions spaces and a map so we can use the Banach Fixed-Point Theorem. In the proceeding sections, we will always assume 0 < r < 1. LetH∞(B(0, r)2) represent the space of all holomorphic functions of two complex variables bounded on ∂B(0, r)2 with the usual norm ‖ · ‖H∞ := sup (z1,z2)∈B(0,r)2 |(·)(z1, z2)|. (3.4) Definition 3.2. Define the subset Hr of H∞(B(0, r)2) such that for each g ∈ Hr, (1) ck,j(g) ≥ 0 for integers k, j ≥ 0 and (2) ∞∑ k=0 ∞∑ j=0 ck,j(g) ≤ 1, (3.5) where ck,j(g) = 1 k!j! ∂k+jg ∂xk∂yj (0, 0) = 1 (2πi)2 ∫ z1∈C(0,r) ∫ z2∈C(0,r) g(z1, z2) zk+1 1 zj+1 2 dz2 dz1, (3.6) is the coefficient of xkyj of g in its power series expansion centered at 0, where x and y are the function parameters. Here C(0, r) represents the path in C along a circle centered at the origin of radius r traversed once counter-clockwise. Proposition 3.3. Hr is a separable, complete metric space with respect to the norm ‖ · ‖H∞ . Proof. To show completeness, we need prove is that Hr is closed in H∞(B(0, r)2). By equation (3.6) it follows that the functions cj,k are continuous on H∞, and hence 6 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 the only difficulty is to show Definition 3.2, Part (3.5) is preserved by convergent sequences in H∞. But this follows by noting that it is equivalent to M∑ k,j=0 ck,j(g) ≤ 1 for all M ∈ N (3.7) Separability follows from part (3.5) since it follows that the set of polynomials in Hr is dense in Hr. � Definition 3.4. Define P[0,∞) to be the set of all partitions of [0,∞). We say a function f : [0,∞)→ H∞(B(0, r)2) has bounded variation on [0,∞) if V (f) := sup S∈P[0,∞) M∑ k=1 ‖f(tk)− f(tk−1)‖H∞ (3.8) is finite, where S is a partition 0 ≤ t0 ≤ t1 ≤ · · · ≤ tM <∞. For a function f : [0,∞) → H∞(B(0, r)2) we define the usual norm using the essential supremum, ‖f‖L∞(H∞) := ess sup t∈[0,∞) ‖f(t)‖H∞ . (3.9) Definition 3.5. Define BV ([0,∞), H∞(B(0, r)2)) to be the space of all functions mapping from [0,∞) to H∞(B(0, r)2) with bounded variation on [0,∞) with norm ‖ · ‖BV (H∞) := ‖ · ‖L∞(H∞) + V (·). (3.10) The following result can be proved using standard techniques (see [9]). Proposition 3.6. BV ([0,∞), H∞(B(0, r)2)) is a separable Banach space, consist- ing of bounded Borel measurable functions from [0,∞) to H∞(B(0, r)2). Thus for any distribution P , we can define integration for a measurable function f : [0,∞) → H∞(B(0, r)2) using the Bochner integral on the space L1(([0,∞), dP ), H∞(B(0, r)2)) [2]. Definition 3.7. Define Br,m to be the subspace of BV ([0,∞), H∞(B(0, r)2)) such that for all f ∈ Br,m, (1) f(t) ∈ Hr for all t ∈ [0,∞), (2) c0,0(f(t)) = 0 for all t ∈ [0,∞), and (3) V (f) ≤ m. Proposition 3.8. Br,m is a complete metric space with the norm ‖ · ‖BV (H∞), and thus the product space (Br,m) 2 is a complete metric space with product norm ‖ ((·)1, (·)2) ‖BV 2(H∞) := ‖(·)1‖BV (H∞) + ‖(·)2‖BV (H∞). (3.11) Definition 3.9. Define T : (Br,m) 2 → (BV ([0,∞), H∞(B(0, r)))) 2 (3.12) to be a map such that for (f, g) ∈ (Br,m)2, (T (f, g))(t, s) := ( (·)1(1− P (t)) + ∫ t 0 f(t− τ)g(t− τ) dP (τ), (·)2(1−Q(s)) + ∫ s 0 f(s− τ)g(s− τ) dQ(τ) ) (3.13) EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 7 for all (t, s) ∈ [0,∞)2, where (·)1 and (·)2 represent the first and second parameters of functions in H∞(B(0, r)2), respectively. Theorem 3.10 (Banach Fixed-Point Theorem). Suppose X is a complete metric space with distance function d(·, ·) and T : X → X is a map such that there exists a constant 0 ≤ γ < 1 where d(T (x), T (y)) ≤ γd(x, y) for all x, y ∈ X. Then T has a unique fixed-point. Theorem 3.11. If 0 < m < 1 2 and 0 < r < min ( 1−2m 4 , m 2(1+m) ) , then T is a contraction mapping of ( Br,m )2 to itself, and hence T has a unique fixed point. Remark 3.12. Since analytic functions are uniquely determined, as long as r > 0, we can extend our fixed point to converge on B(0, 1)2. Theorem 3.13. Let (f, g) be the fixed point of T . Then f(t, 1, 1) = g(t, 1, 1) = 1, and hence f and g are generating functions that solve equation (3.3). Proof. Clearly f(t, 1, 1), g(t, 1, 1) ≤ 1. Next, since P (0) = Q(0) = 0, we see that f(0, 1, 1) = g(0, 1, 1) = 1. (3.14) Let t∗ = inf{t ≥ 0 : f(t, 1, 1) 6= 1 or g(t, 1, 1) 6= 1}. (3.15) Pick 0 < ε < 1. There exists δ > 0 such that P (δ), Q(δ) < ε. Let f0 = inf{f(t, 1, 1), g(t, 1, 1) : t ∈ [t∗, t∗ + δ)}. (3.16) From equation (3.3), we obtain that for t ∈ [t∗, t∗ + δ) f(t, 1, 1) = 1− P (t) + ∫ t−t∗<τ≤t f(t− τ, 1, 1)g(t− τ, 1, 1) dP (τ) + ∫ t−t∗ 0 f(t− τ, 1, 1)g(t− τ, 1, 1) dP (τ) ≥ 1− P (t− t∗) + P (t− t∗)f20 (3.17) and similarly g(t, 1, 1) ≥ 1−Q(t− t∗) +Q(t− t∗)f20 . (3.18) Hence f0 ≥ 1− ε+ εf20 ⇒ (f0 − 1)(εf0 − 1 + ε) ≤ 0 (3.19) from which it follows that f0 ≥ 1. � 3.4. Series solutions to the generating function equation. If we are only con- cerned about the total cell population and not specifically the mother and daughter cell populations, we can write our generating functions f(t, x, y) and g(t, x, y) as simply f(t, x) = f(t, x, x) and g(t, x) = g(t, x, x). So, we can express our generating functions f(t, x) and g(t, x) as series f(t, x) = ∞∑ k=0 ck(t)xk, g(t, x) = ∞∑ k=0 bk(t)xk, (3.20) which are necessarily convergent for all t ≥ 0 when x ∈ B(0, 1). If f(t, x) and g(t, x) satisfy our integral equations (3.3), then by matching up coefficients and 8 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 noting c0(t) = b0(t) = 0 for all t ≥ 0, we obtain ck(t) = { 1− P (t) k = 1∫ t 0 ∑k−1 j=1 cj(t− τ)bk−j(t− τ) dP (τ) k ≥ 2 bk(t) = { 1−Q(t) k = 1∫ t 0 ∑k−1 j=1 cj(t− τ)bk−j(t− τ) dQ(τ) k ≥ 2 . (3.21) In Section 6, we will use these formulas to determine the distribution of the mutant cell population when P and Q are multi-phase distributions. 4. Asymptotics of cell growth Definition 4.1. A probability distribution P is a δ-lattice distribution if P is constant except at jumps at multiples of some δ > 0. Define the following integrals for s ∈ C: p∗(s) := ∫ ∞ 0 e−st dP (t), q∗(s) := ∫ ∞ 0 e−st dQ(t), ψ(s) := 1− p∗(s)− q∗(s). Note that ψ′(s) = ∫ ∞ 0 te−st d(P (t) +Q(t)), (4.1) is positive when s ≥ 0 is real. Proposition 4.2. ψ has a unique real root, α ∈ (0,∞). Proof. This follows since ψ(0) = −1, lims→+∞ ψ(s) = 1, and ψ′(s) > 0 for s in [0,∞). � Remark 4.3. In the rest of this paper, unless otherwise stated, α will always refer to the unique real root of ψ. The goal of Section 4 is to prove the following result. Theorem 4.4. If for any δ > 0 we have that either P or Q is not a δ-lattice distribution, then there exists a non-negative random variable V such that Xte −αt → p∗(α)V, Yte −αt → q∗(α)V (4.2) where the convergence is in L2. 4.1. The expectations E[Xt] and E[Yt]. Define the following expectations: mf (t) := E [Xt | X0 = 1, Y0 = 0] = fx(t, 1, 1); mg(t) := E [Xt | X0 = 0, Y0 = 1] = gx(t, 1, 1); nf (t) := E [Yt | X0 = 1, Y0 = 0] = fy(t, 1, 1); ng(t) := E [Yt | X0 = 0, Y0 = 1] = gy(t, 1, 1). (4.3) EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 9 Proposition 4.5. Expectations mf (t), mg(t), nf (t), and ng(t) satisfy the system mf (t) = 1− P (t) + ∫ t 0 mf (t− τ) +mg(t− τ) dP (τ) mg(t) = ∫ t 0 mf (t− τ) +mg(t− τ) dQ(τ) nf (t) = ∫ t 0 nf (t− τ) + ng(t− τ) dP (τ) ng(t) = 1−Q(t) + ∫ t 0 nf (t− τ) + ng(t− τ) dQ(τ). (4.4) Proof. Suppose x, y ∈ B(0, 1). We can differentiate with respect to x and y, both sides of both equations in (3.3). It is straightforward to bring the derivatives under the integral sign using the Cauchy integral formula and Fubini’s Theorem. Since these generating functions and their derivatives have positive coefficients and are increasing in both the x and y parameters, we can let x, y → 1−. � Proposition 4.6. Expectations mf (t), mg(t), nf (t), and ng(t) are non-decreasing and finite for t ≥ 0. Proof. Since no cells can die, mf , mg, nf , and ng must be non-decreasing. Since P and Q are life-span distributions with P (0) = Q(0) = 0 and P and Q are right continuous, then for ε > 0, there exists a δ > 0 such that P (δ) < ε and Q(δ) < ε. Now, rewriting and bounding system (4.4), we obtain mf (t) ≤ 1 + εmf (t) +mf (t− δ) + εmg(t) +mg(t− δ) mg(t) ≤ εmf (t) +mf (t− δ) + εmg(t) +mg(t− δ). (4.5) So, mf (t) +mg(t) ≤ 1 + 2ε(mf (t) +mg(t)) + 2(mf (t− δ) +mg(t− δ)), (4.6) and mf (t) +mg(t) ≤ 1 + 2(mf (t− δ) +mg(t− δ)) 1− 2ε , (4.7) that is, if mf (t− δ) and mg(t− δ) are finite, then mf (t) and mg(t) are finite. Since mf (0) +mg(0) = 1, we can inductively conclude mf (t) and mg(t) are finite for all t ≥ 0, and the result is proven for mf and mg. nf and ng are similar. � 4.2. Convergence of E[Xt]e −αt and E[Yt]e −αt. We will compute Laplace trans- forms, and then analyze their poles. We need to first review some properties of the Laplace transform. Definition 4.7. Let h : [0,∞)→ R. The Laplace transform of h is (Lh)(s) := ∫ ∞ 0 e−sth(t) dt, s ∈ C. (4.8) We will often write h∗ := Lh to represent the Laplace transform of h. Proposition 4.8 (see [13, page 92]). If f, g ∈ L1([0, R)) for all R > 0, then L(f ∗ g) := (Lf)(Lg), (4.9) provided all three transforms exist. Similarly we also have the following result. 10 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 Proposition 4.9. Let ϕ : [0,∞) → R be a function and G be a probability distri- bution on [0,∞). Define (ϕ ∗ dG)(t) := ∫ t 0 ϕ(t− τ) dG(τ), (4.10) g∗(s) := ∫ ∞ 0 e−st dG(t). (4.11) Then L(ϕ ∗ dG) = (Lϕ)g∗ (4.12) when all integrals converge. Let s ∈ C with Re(s) > 0. Then m∗f (s) := (Lmf )(s) = 1 s − 1 s p∗(s) +m∗f (s)p∗(s) +m∗g(s)p ∗(s), m∗g(s) := (Lmg)(s) = m∗f (s)q∗(s) +m∗g(s)q ∗(s), (4.13) Solving for m∗f (s) and m∗g(s), we obtain m∗f (s) = 1− p∗(s)− q∗(s) + p∗(s)q∗(s) s(1− p∗(s)− q∗(s)) , m∗g(s) = (1− p∗(s))q∗(s) s(1− p∗(s)− q∗(s)) . (4.14) Similarly, n∗f (s) = (1− q∗(s))p∗(s) s(1− p∗(s)− q∗(s)) , n∗g(s) = 1− p∗(s)− q∗(s) + p∗(s)q∗(s) s(1− p∗(s)− q∗(s)) . (4.15) If P and Q are lattice distributions, then it can be shown that there are infinitely many zeros of ψ with Re(s) = α, but if P or Q is not a lattice distribution, then the following hold. Proposition 4.10. If P and Q are not both δ-lattice distributions, and α + iτ is a zero of ψ, then τ = 0. Hence m∗f , m∗g, n∗f , and n∗g are analytic at α + iτ when τ 6= 0. Proof. Since α and α+ iτ are zeros, we have Re (ψ(α)− ψ(α+ iτ)) = ∫ ∞ 0 e−αt(cos(τt)− 1) d(P (t) +Q(t)) = 0, (4.16) So, ∫ ∞ 0 e−αt(cos(τt)− 1) d(P (t) +Q(t)) = 0, (4.17) and since e−αt(cos(τt)− 1) ≤ 0 (4.18) for all t ≥ 0, we can conclude that cos(τt)− 1 = 0, P +Q-a.e., (4.19) which requires τ = 0 or t ∈ 2π τ Z. But if it’s the latter, then P +Q, and hence P and Q, are necessarily 2π τ -lattice distributions, which is a contradiction. So τ = 0. � EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 11 Define the following constants: c1 = p∗(α)q∗(α) αψ′(α) , c2 = [q∗(α)]2 αψ′(α) , d1 = [p∗(α)]2 αψ′(α) . (4.20) Proposition 4.11. m∗f , m∗g, n∗f , and n∗g have a poles at α of order 1 with residues c1, c2, d1, c1 respectively. Now, we use the Wiener-Ikehara Tauberian Proposition to show m converges to an exponential function if P is not a lattice distribution. Theorem 4.12 (Wiener-Ikehara Theorem [13, p. 233]). If ϕ(t) is a non-negative, non-decreasing function for t ≥ 0 such that the integral f(s) = ∫ ∞ 0 e−stϕ(t) dt, s = σ + iτ ∈ C (4.21) converges for σ > 1, and if for some constants A ∈ C, α > 0, and some function g(τ), lim σ→1+ f(s)− A s− α = g(τ) (4.22) uniformly in every finite interval −a ≤ τ ≤ a, then lim t→∞ ϕ(t)e−αt = A. (4.23) Note that the result is stated for α = 1 in [13], but the general case is easily seen to follow. Theorem 4.13. Suppose P and Q are not both δ-lattice distributions. Then mf (t) ∼ c1eαt, mg(t) ∼ c2eαt, nf (t) ∼ d1eαt, ng(t) ∼ c1eαt. (4.24) 4.3. Convergence of Xt/E[Xt] and Yt/E[Yt]. Let F represent the joint generating function of Xt, Yt, Xt+τ , and Yt+τ when you start with exactly one mother cell and no daughter cells, and let G represent the joint generating function when you start with exactly one daughter cell and no mother cells. Then F (t, τ, x1, x2, y1, y2) = x1x2(1− P (t+ τ)) + ∫ t 0 F (t− y, τ, x1, x2, y1, y2)G(t− y, τ, x1, x2, y1, y2) dP (y) + x1 ∫ τ t f(t+ τ − y, x2, y2)g(t+ τ − y, x2, y2) dP (y) G(t, τ, x1, x2, y1, y2) = y1y2(1−Q(t+ τ)) + ∫ t 0 F (t− y, τ, x1, x2, y1, y2)G(t− y, τ, x1, x2, y1, y2) dQ(y) y1 ∫ τ t f(t+ τ − y, x2, y2)g(t+ τ − y, x2, y2) dQ(y), where the parameters x1 and y1 correspond to to Xt and Yt, and x2 and y2 corre- spond to Xt+τ and Yt+τ . 12 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 Define the following expectations: m2,f (t, τ) := E[XtXt+τ | X0 = 1, Y0 = 0] = Fx1x2 (t, τ, 1, 1, 1, 1); m2,g(t, τ) := E[XtXt+τ | X0 = 0, Y0 = 1] = Gx1x2 (t, τ, 1, 1, 1, 1); n2,f (t, τ) := E[YtYt+τ | X0 = 1, Y0 = 0] = Fy1y2(t, τ, 1, 1, 1, 1); n2,g(t, τ) := E[YtYt+τ | X0 = 0, Y0 = 1] = Gy1y2(t, τ, 1, 1, 1, 1); cf (t) := E[XtYt | X0 = 1, Y0 = 0] = Fx1y1(t, t, 1, 1, 1, 1); cg(t) := E[XtYt | X0 = 0, Y0 = 1] = Gx1y1(t, t, 1, 1, 1, 1). (4.25) Proposition 4.14. The expectations m2,f and m2,g satisfy the system m2,f (t, τ) = 1− P (t+ τ) + ∫ t 0 m2,f (t− y, τ) +m2,g(t− y, τ) dP (y) + ∫ t 0 mf (t− y)mg(t+ τ − y) +mf (t+ τ − y)mg(t− y) dP (y) + ∫ t+τ t mf (t+ τ − y) +mg(t+ τ − y) dP (y), (4.26) m2,g(t, τ) = ∫ t 0 m2,f (t− y, τ) +m2,g(t− y, τ) dQ(y) + ∫ t 0 mf (t− y)mg(t+ τ − y) +mf (t+ τ − y)mg(t− y) dQ(y). (4.27) Similarly, the expectations n2,f and n2,g satisfy the system n2,f (t, τ) = ∫ t 0 n2,f (t− y, τ) + n2,g(t− y, τ) dP (y) + ∫ t 0 nf (t− y)ng(t+ τ − y) + nf (t+ τ − y)ng(t− y) dP (y), (4.28) n2,g(t, τ) = 1−Q(t+ τ) + ∫ t 0 n2,f (t− y, τ) + n2,g(t− y, τ) dQ(y) + ∫ t 0 nf (t− y)mg(t+ τ − y) + nf (t+ τ − y)mg(t− y) dQ(y) + ∫ t+τ t nf (t+ τ − y) + ng(t+ τ − y) dQ(y), and the expectations cf and cg satisfy the system cf (t) = ∫ t 0 cf (t− y) + cg(t− y) dP (y) + ∫ t 0 mf (t− y)ng(t− y) + nf (t− y)mg(t− y) dP (y), cg(t) = ∫ t 0 cf (t− y) + cg(t− y) dQ(y) + ∫ t 0 mf (t− y)ng(t− y) + nf (t− y)mg(t− y) dQ(y). EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 13 Following the arguments given in Proposition 4.6, we obtain the following result. Proposition 4.15. The expectations m2,f (t, τ), m2,g(t, τ), n2,f (t, τ), n2,g(t, τ) are non-decreasing in both arguments and finite for t, τ ≥ 0. Similarly cf (t) and cg(t) are non-decreasing and finite for t ≥ 0. Lemma 4.16 ([1, Lemma 2]). If v(t) satisfies the equation v(t) = ∫ t 0 v(t− y)dH(y) + h(t) (4.29) where H is a non-decreasing function with H(0) = 0 and H(∞) = α < 1, and h(t) is a bounded function such that limt→∞h(t) = c, then lim t→∞ v(t) = c 1− α (4.30) Now, we can prove the following convergence theorem about E[XtXt+τ ] and E[YtYt+τ ]. Theorem 4.17. If P and Q is not both δ-lattice distributions, then m2,f (t, τ) ∼ D1e αteα(t+τ), n2,f (t, τ) ∼ D2e αteα(t+τ), (4.31) uniformly in τ , and cf (t) ∼ D1e 2αt, (4.32) where D1 = c1c2p ∗(2α) ψ(2α) , D2 = c1c2q ∗(2α) ψ(2α) . (4.33) There are similar results for m2,g, n2,g and cg. Proof. Multiply both sides of (4.26) by e−αte−α(t+τ), and multiply both sides of (4.27) by e−αte−α(t+τ), and set K1(t, τ) := e−αte−α(t+τ)m2,f (t, τ), K2(t, τ) := e−αte−α(t+τ)m2,g(t, τ), (4.34) dP (y) := e−2αydP (y) p∗(2α) , dQ(y) := e−2αydQ(y) q∗(2α) (4.35) to obtain K1(t, τ) = p∗(2α) ∫ t 0 K1(t− y, τ) +K2(t− y) dP (y) + h1(t, τ), K2(t, τ) = q∗(2α) ∫ t 0 K1(t− y, τ) +K2(t− y) dQ(y) + h2(t, τ) (4.36) where h1(t, τ) :=e−αte−α(t+τ)(1 + P (t+ τ)) + p∗(2α) ∫ t 0 e−α(t−y)mf (t− y)e−α(t+τ−y)mg(t+ τ − y) dP (y) + p∗(2α) ∫ t 0 e−α(t−y)mg(t− y)e−α(t+τ−y)mf (t+ τ − y) dP (y) + p∗(2α) ∫ t+τ t e−α(t−y)e−α(t+τ−y)mf (t+ τ − y) dP (y) + p∗(2α) ∫ t+τ t e−α(t−y)e−α(t+τ−y)mg(t+ τ − y) dP (y) (4.37) 14 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 and h2(t, τ) :=q∗(2α) ∫ t 0 e−α(t−y)mf (t− y)e−α(t+τ−y)mg(t+ τ − y) dQ(y) + q∗(2α) ∫ t 0 e−α(t−y)mg(t− y)e−α(t+τ−y)mf (t+ τ − y) dQ(y). Using Theorem 4.13, a standard ε-δ argument shows that lim t→∞ h1(t, τ) = 2p∗(2α)c1c2, lim t→∞ h2(t, τ) = 2q∗(2α)c1c2 (4.38) uniformly in τ . Now, adding Equations (4.36) and setting K(t, τ) := K1(t, τ) +K2(t, τ), h(t, τ) := h1(t, τ) + h2(t, τ), dR(y) := e−2αyd(P (y) +Q(y)) p∗(2α) + q∗(2α) , (4.39) we obtain K(t, τ) = (p∗(2α) + q∗(2α)) ∫ t 0 K(t− y, τ) dR(y) + h(t, τ). (4.40) Since p∗ + q∗ is a decreasing function and p∗(α) + q∗(α) = 1, then we must have p∗(2α) + q∗(2α) < 1. Moreover, from (4.38), lim t→∞ h(t, τ) = 2(p∗(2α) + q∗(2α))c1c2 (4.41) uniformly in τ . Using Lemma 4.16, we can conclude we can conclude that lim t→∞ K(t, τ) = 2(p∗(2α) + q∗(2α))c1c2 ψ(2α) . (4.42) Looking back at (4.36), and using (4.38) and (4.42), we obtain the equation for mf (t, τ). The other equations follow similarly. � Now, we define Wt := Xt c1eαt , Vt := Yt d1eαt . (4.43) Proposition 4.18. If P and Q are not both δ-lattice distributions, then lim t→∞ E [ (Wt+τ −Wt) 2 ] = 0, lim t→∞ E [ (Vt+τ − Vt)2 ] = 0 (4.44) uniformly in τ , and E[(Wt − Vt)2]→ 0 (4.45) Proof. First we show these results conditionally on either X0 = 1 and Y0 = 0, or X0 = 0 and Y0 = 1. The case when X0 = 1 and Y0 = 0 follows by multiplying out E [ (Wt+τ −Wt) 2 ] and applying Theorem 4.17 to mf (t, t), mf (t, τ), and mf (τ, τ). The other cases follow similarly. Now, suppose that X0 and Y0 are not specified. Then we can create X0 random variables (X̃ (n) t , Ỹ (n) t ) for 1 ≤ n ≤ X0, each being an independent copy of (Xt, Yt) conditioned on X0 = 1, Y0 = 0, and Y0 random variables (X̂ (n) t , Ŷ (n) t ) for 1 ≤ n ≤ EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 15 Y0, each being an independent copy of (Xt, Yt) conditioned on X0 = 0, Y0 = 1. Then Xt = X0∑ n=1 X̃ (n) t + Y0∑ n=1 X̂ (n) t Yt = X0∑ n=1 Ỹ (n) t + Y0∑ n=1 Ŷ (n) t (4.46) A short verification shows that the result still holds. � Since L2 is a complete metric space, we obtain the immediate corollary. Corollary 4.19. Wt and Vt converge in L2 to a random variable W . 5. Distribution of mutant cells In this section, we will determine a probability distribution for the mutant cell population size when the total cell population is effectively infinite. We start by making assumptions about the cell population and the expected number of muta- tions that occur at any given time. The culmination of this paper are Theorems 5.8 and 5.9, where we give an explicit formula for the generating functions for the mutant cell population. 5.1. Preliminaries and assumptions. We start by stating our assumptions about cell growth under asymmetric cell division in which cells can mutate. Assumption 5.1. Cells have the following properties: (1) there are exactly two types of cells: mother cells and daughter cells; (2) all cells are independent of each other, mother cells are identical to other mother cells, and daughter cells are identical to other daughter cells; (3) cell life-span for mother cells and daughter cells are strictly positive, real- valued random variables Tx and Ty, respectively, with distributions P and Q, respectively; (4) at the end of a cell’s life, both mother cells and daughter cells will divide into one mother cell and one daughter cell; (5) when a cell divides, exactly one child cell can mutate with probability µ; (6) all children spawned from a mutant cell will be mutants with no chance of losing the mutation. We will be working backwards in time and let the current (time 0) cell population be n. In order for our model to work, we have to assume the following about our mutation rate and current cell population: Assumption 5.2. We will assume n is very large and µ is sufficiently small so the product of µ and n stays fixed as µ→ 0+ and n→∞, with m := µn (5.1) 5.2. Constructing the generating function. We will start by assuming time is discrete and takes values kδt for integers k ≥ 0, where δt > 0. Because we are working backwards in time, kδt represents time k time-units ago. Consider 4 random variables: (1) the number of mutant mother cells NM k created k time-units ago; (2) the number of mutant daughter cells ND k created k time-units ago; 16 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 (3) the number of mutants Mk that arise from any one mutant mother cell which was created at k time-units ago; (4) the number of mutants Dk that arise from any one mutant daughter cell which was created at k time-units ago. Then the total number of mutants will be given by the formula R = RM +RD, (5.2) where RM = ∞∑ k=0 NM k∑ j=1 M (j) k , RD = ∞∑ k=0 ND k∑ j=1 D (j) k , (5.3) and where M (j) k and D (j) k denote independent copies of Mk and Dk, respectively. Let Fm(t) and FD(t) be functions representing the mother and daughter cell population t time-units ago, respectively. So, FM (0) + FD(0) = n. Assumption 5.3. We will assume NM k is a Poisson random variable such that NM k ∼ Pois (µ (FM (kδt)− FM ((k + 1)δt))) , (5.4) which is assuming on average, a µ-proportion of new mother cells created during that time-unit will become mutants. Similarly, we will assume ND k is a Poisson random variable such that ND k ∼ Pois (µ (FD(kδt)− FD((k + 1)δt))) , (5.5) which is assuming on average, a µ-proportion of new daughter cells created during that time-unit will become mutants. Remark 5.4. In practice, k in Equations (5.3) only needs to increase to a point where the population of non-mutants k time-units is much smaller than the current population. However, since we will assume we start with an effectively infinite cell population, mathematically, letting k →∞ is reasonable. Assumption 5.5. In our model, we are counting the cell population that grows from a single mutant cell, and because cell growth is independent, we can assume RM and RD are independent. Using Assumption 5.3, we obtain the generating function for RM , GRM (x) = ∞∏ k=0 GNM k (GMk (x)) = ∞∏ k=0 exp (µ (FM (kδt)− FM ((k + 1)δt)) (GMk (x)− 1)) = exp ( − µFM (0) + µ ∞∑ k=0 (FM (kδt)− FM ((k + 1)δt))GMk (x) ) . (5.6) Letting δt → 0 and replacing GMk (x) with f(t, x), we obtain a continuous-time model, gRM (x) = exp ( − µFM (0)− µ ∫ ∞ 0 f(t, x) dFM (t) ) . (5.7) Similarly, gRD (x) = exp ( − µFD(0)− µ ∫ ∞ 0 g(t, x) dFD(t) ) . (5.8) EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 17 Remark 5.6. f(t, x) and g(t, x) represent the generating function for the total cell population when you start from one mother cell or one daughter cell, respectively. This means f(t, x) = f(t, x, x) and g(t, x) = g(t, x, x) where the right-side of the equalities are the generating functions in (3.3). Using Assumption 5.5, we obtain gR(x) = gRM (x)gRD (x) = exp ( − µn− µ ∫ ∞ 0 f(t, x) dFM (t)− µ ∫ ∞ 0 g(t, x) dFD(t) ) . (5.9) Since we are working backwards in time, we can interpret our current time as infinite relative to when we started from a single cell. Suppose P and Q are not both δ-lattice distributions. By Theorem 4.4, we know that asymptotically, the proportion of mother cells to the total cell population will be p∗(α), while the proportion of daughter cells to the total cell population will be q∗(α) = 1− p∗(α). This leads us to the following assumption about the cell population at backwards time t. Assumption 5.7. The mother cell population at backwards time t is FM (t) = np∗(α)e−αt, (5.10) and the daughter cell population at backwards time t is Fd(t) = nq∗(α)e−αt (5.11) with high probability. Hence gR(x) = exp ( −m+mα ∫ ∞ 0 (p∗(α)f(t, x) + q∗(α)g(t, x))e−αt dt ) , (5.12) where we recall Equation (5.1). The above construction culminates to the following theorem. Theorem 5.8. Suppose distributions P and Q are not both δ-lattice distributions for any δ > 0. Let R be a random variable representing the mutant cell population in an effectively infinite cell population n with effectively zero mutation rate µ, and set m := µn. Let f(t, x) and g(t, x) be the generating functions for the total cell population at time t when you start from a single mother cell and daughter cell, respectively, with life-span distributions P and Q for mother cells and daughter cells, respectively, which must satisfy f(t, x) = x(1− P (t)) + ∫ t 0 f(t− τ, x)g(t− τ, x) dP (τ) g(t, x) = x(1−Q(t)) + ∫ t 0 f(t− τ, x)g(t− τ, x) dQ(τ). (5.13) Let α be the root of 1− ∫ ∞ 0 e−st dP (t)− ∫ ∞ 0 e−st dQ(t) = 0 (5.14) and γ := ∫ ∞ 0 e−αt dP (t). (5.15) 18 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 Then the generating function for R is gR(x) = exp ( −m+mα ∫ ∞ 0 (γf(t, x) + (1− γ)g(t, x))e−αt dt ) . (5.16) It is also worth stating the special case when cell division is symmetric, when we do not need to distinguish between mother and daughter cells. Theorem 5.9. Suppose the distribution P is not a δ-lattice distribution for any δ > 0. Let R be a random variable representing the mutant cell population in an effectively infinite cell population n with effectively zero mutation rate µ, and set m := µn. Let f(t, x) total cell population at time t when you start from a cell with life-span distribution P , which must satisfy f(t, x) = x(1− P (t)) + ∫ t 0 f(t− τ, x)2 dP (τ) (5.17) Let α be the root of 1− 2 ∫ ∞ 0 e−st dP (t) = 0. (5.18) Then the generating function for R is gR(x) = exp ( −m+mα ∫ ∞ 0 f(t, x)e−αt dt ) . (5.19) 6. Examples 6.1. Life-span is multi-phase. We say that P has a “multi-phase” distribution with parameter n if is a gamma distribution with parameters β = 1 and α = n. Biologically, this means that a cell is more likely to divide around the nth “stage” of its life. If n = 1, then a cell is most likely to divide closer to birth, and P would be an exponential distribution. However, when n = 3, for example, it would mean the cell takes time to mature before it can divide, which is more realistic biologically. Suppose P and Q are multi-phase distributions with parameter n and m, respec- tively. Then dP (t) = tk−1e−t (k − 1)! dt, dQ(t) = tn−1e−t (n− 1)! dt (6.1) and the exponential growth rate α for both Xt and Yt is the root to (1 + s) k+n − (1 + s) k − (1 + s) n = 0. (6.2) The multi-phase distribution is of great importance in the asymmetric cell di- vision case, and is a major inspiration for this paper, where we were interested in developing a generating function when the daughter cells take time to mature before they can divide. A good model for this would be to assume mother cells have multi-phase distribution with parameter k = 1 or n = 2, while daughter cells have multi-phase distribution k = 2 or n = 3. Theorem 6.1. If P and Q are multi-phase distributions with parameters k and n, respectively, then the generating functions f = f(t, x, y) and g = g(t, x, y) satisfy( 1 + ∂ ∂t )k f = fg, ( 1 + ∂ ∂t )n g = fg, f(0, x, y) = x, g(0, x, y) = y (6.3) EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 19 6.2. Symmetric cell division. Suppose P is an exponential distribution with pa- rameter 1, that is, it is a multiphase distributions with n = 1. Then the differential equation for the generating function is ∂f ∂t = f2 − f, f(0, x) = x (6.4) which can be solved to obtain f(t, x) = xe−t 1− x+ xe−t . (6.5) Moreover α = 1. So, the generating function becomes gR(x) = (1− x) m(1−x) x = e−m + 1 2 e−mmx+ 1 24 e−mm(3m+ 4)x2 + 1 48 e−mm(m+ 2)2x3 + . . . . which is the Lea-Coulson distribution (see [5]). 6.3. Asymmetric cell division. Suppose P and Q are multi-phase distributions with parameters 1 and 2, respectively. Then we have α = γ = √ 5− 1 2 . (6.6) This models a situation where mother cells have shorter expected life-spans than daughter cells. Using the series solution to f(t, x) and g(t, x) discussed in Section 3.4, and the Ma-Sandri-Sarkar algorithm (see [7]), we can compute the coefficients of gR(x) relatively quickly. These were computed using the software package sagemath (see [10]), using code shown in Section 6.4. Pr(R = 0) = e−m Pr(R = 1) = 0.472135955me−m Pr(R = 2) = (0.111456180002m2 + 0.172209268743m)e−m Pr(R = 3) = (0.0175408233286m3 + 0.0813061875578m2 + 0.0923265707329m)e−m Pr(R = 4) = (0.00207041334343m4 + 0.019193787255m3 + 0.0584187097653m2 + 0.054750171345m)e−m Pr(R = 5) = (0.000195503316229m5 + 0.00302069235857m4 + 0.0172912064384m3 + 0.0417490156659m2 + 0.0359263840576m)e−m Pr(R = 6) = (1.53840241522× 10−5m6 + 0.000356544367868m5 + 0.00327215809956m4 + 0.0144601403668m3 + 0.0306527224493m2 + 0.0253788856813m)e−m Pr(R = 7) = (1.03762156212× 10−6m7 + 3.36674831246× 10−5m6 + 0.000451249775082m5 + 0.00313432945934m4 + 0.0118370547378m3 + 0.0232240363224m2 20 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 + 0.0188813484058m)e−m Pr(R = 8) = (6.12373058949× 10−8m8 + 2.64927154958× 10−6m7 + 4.87502731868× 10−5m6 + 0.000487114247624m5 + 0.00283908715808m4 + 0.00968207809151m3 + 0.018100793276m2 + 0.0145941835643m)e−m Pr(R = 9) = (3.21248154448× 10−9m9 + 1.78688050445× 10−7m8 + 4.3192805187× 10−6m7 + 5.85526505241× 10−5m6 + 0.000484413747255m5 + 0.00250654733417m4 + 0.0079759775165m3 + 0.0144521031634m2 + 0.0116173054839m)e−m Pr(R = 10) = (1.51672804192× 10−10m10 + 1.0545631668× 10−8m9 + 3.23914598676× 10−7m8 + 5.71296309855× 10−6m7 + 6.34669453811× 10−5m6 + 0.000460389737452m5 + 0.00218929525915m4 + 0.00663573155793m3 + 0.0117763023231m2 + 0.00946644120115m)e−m 6.4. Sage code to compute probabilities for multiphase distributions. t = var(’t’); m = var(’m’); y = var(’y’); a = var(’a’) # The next line seems to be required to make the definite integrals # work. assume(t>0) # number of coefficients to compute. degree = 10 # life-span distribution density of P (multi-phase with parameter # k). k = 1 p(t) = t^(k-1)*exp(-t)/factorial(k-1) P(t) = integrate(p(y),y,0,t) # life-span distribution density of Q (multi-phase with parameter # n). n = 2 q(t) = t^(n-1)*exp(-t)/factorial(n-1) Q(t) = integrate(q(y),y,0,t) # find alpha and gamma. eqn = (1+x)^(k+n)-(1+x)^k-(1+x)^n assume(x,’real’) alpha = max(a.rhs() for a in eqn.solve(x)) # use float for faster calculations. EJDE-2021/56 BRANCHING PROCESSES AND THE LURIA-DELBRÜCK EXPERIMENT 21 alpha = float(alpha) gamma = integrate(exp(-alpha*t)*p(t),t,0,infinity) # store coefficients c[i] and b[i] of the generating functions f and # g. c = [0,1-P(t)] b = [0,1-Q(t)] for i in xrange(2,degree+1): tempsum(t) = sum(expand(c[j]*b[i-j]) for j in xrange(1,i)) # indefinite integration is faster than definite integration. temp1(y)=integrate(expand(tempsum(y)*expand(p(t-y))),y) temp2(y)=integrate(expand(tempsum(y)*expand(q(t-y))),y) c.append(temp1(t)-temp1(0)) b.append(temp2(t)-temp2(0)) # store coefficients of # H = integral(pgf*alpha*e^(-alpha t),t,0,infinity). h = [0] for i in xrange(1,degree+1): h.append(integrate((gamma*c[i]+ (1-gamma)*b[i])* alpha*exp(-alpha*t),t,0,infinity)) # compute coefficients of e^H using MSS algorithm. prob = [1] for r in xrange(1,degree+1): prob.append(expand(sum(m*s/r*h[s]*prob[r-s] for s in xrange (1,r+1)))) for r in xrange(0,degree+1): print "\\Pr(R =",r,") &=",latex(e^(-m)*prob[r]),"\\\\" 6.5. Mathematica code to compute probabilities for multiphase distribu- tions. degree = 10; k = 1; p[t_] = t^(k - 1) Exp[-t]/(k - 1)!; P[t_] = Integrate[p[y], {y, 0, t}]; n = 2; q[t_] = t^(n - 1) Exp[-t]/(n - 1)!; Q[t_] = Integrate[q[y], {y, 0, t}]; alpha = N[Max[x /. Solve[(1 + x)^(k + n) - (1 + x)^k - (1 + x)^n == 0, x]]]; gamma = Integrate[Exp[-alpha t] p[t], {t, 0, Infinity}]; c[1] = 1 - P[t]; b[1] = 1 - Q[t]; Do[tempsum[t_] = Sum[Expand[c[j] b[i - j]], {j, 1, i - 1}]; c[i] = Integrate[Expand[tempsum[y] Expand[p[t - y]]], {y, 0, t}]; b[i] = Integrate[Expand[tempsum[y] Expand[q[t - y]]], {y, 0, t}], 22 S. J. MONTGOMERY-SMITH, H. OVEYS EJDE-2021/56 {i, 2, degree}]; Do[h[i] = Integrate[(gamma c[i] + (1 - gamma) b[i]) alpha Exp[-alpha t], {t, 0, Infinity}], {i, 1, degree}]; prob[0] = 1; Do[prob[r] = Expand[Sum[m s/r h[s] prob[r - s], {s, 1, r}]], {r, 1, degree}]; Table[{r, prob[r] Exp[-m]}, {r, 0, degree}] // TableForm Acknowledgments. We would like to express our sincere thanks to Professor George Smith for bringing this problem to our attention, and many valuable dis- cussions. Some of this research was funded by NSF DMS Grant 0928053, PRISM: Mathematics in Life Sciences. References [1] R. Bellman, T. E. Harris; On age-dependent binary branching processes, Annals of Mathe- matics, 55:280–295, 1952. [2] J. Diestel; Sequences and Series in Banach Spaces, volume 92 of Graduate Texts in Mathe- matics, Springer, 1984. [3] G.B. Folland; Real Analysis: Modern Techniques and Their Applications, Pure and Applied Mathematics: A Wiley Series of Texts, Monographs and Tracts. Wiley, 2013. [4] T. E. Harris; The Theory of Branching Processes, Courier Corporation, 2002. [5] D. E. Lea, C. A. Coulson; The distribution of the numbers of mutants in bacterial populations, Journal of Genetics, 49(3):264–285, 1949. [6] S. E. Luria, M. Delbrück; Mutations of bacteria from virus sensitivity to virus resistance, Genetics, 28(6):491, 1943. [7] W. T. Ma, G. Vh. Sandri, S. Sarkar; Analysis of the Luria-Delbrück distribution using discrete convolution powers, Journal of Applied Probability, pages 255–267, 1992. [8] S. J. Montgomery-Smith, A. Le, G. P. Smith, H. Oveys, D. Pisechko, A. Yates; Estimation of mutation rates from fluctuation experiments via probability generating functions, Unpub- lished, 2015. arXiv:1608.04175. [9] H. Oveys; Age-dependent Branching Processes and Applications to the Luria-Delbrück Ex- periment, PhD thesis, University of Missouri, 2015. [10] The Sage Developers; SageMath, the Sage Mathematics Software System (Version 7.3), 2016, http://www.sagemath.org. [11] S. Sarkar; Haldane’s solution of the Luria-Delbruck distribution, Genetics, 127(2):257, 1991. [12] G. P. Smith, M. Golomb, S. K. Billstein, S. J. Montgomery-Smith; An enduring legacy: The Luria-Delbrück fluctuation test as a classroom investigation in Darwinian evolution, American Biology Teacher, 77(8):614–619, 2015. [13] D. V. Widder; The Laplace Transform, Princeton University Press, 1946. [14] B. Ycart; Fluctuation analysis: can estimates be trusted? PloS one, 8(12):e80958, 2013. Stephen J. Montgomery-Smith Department of Mathematics, University of Missouri, Columbia MO 65211, USA Email address: stephen@missouri.edu Hesam Oveys Courant Institute of Mathematical Sciences, New York University, New York, NY 10012, USA Email address: hoveys@cims.nyu.edu 1. Introduction 1.1. Luria-Delbrück experiment 1.2. The problem 2. Probability generating functions 3. Cell growth under asymmetric division 3.1. Preliminaries and assumptions 3.2. Constructing the generating function 3.3. Existence and uniqueness 3.4. Series solutions to the generating function equation 4. Asymptotics of cell growth 4.1. The expectations E[Xt] and E[Yt] 4.2. Convergence of E[Xt]e-t and E[Yt]e-t 4.3. Convergence of Xt/E[Xt] and Yt/E[Yt] 5. Distribution of mutant cells 5.1. Preliminaries and assumptions 5.2. Constructing the generating function 6. Examples 6.1. Life-span is multi-phase 6.2. Symmetric cell division 6.3. Asymmetric cell division 6.4. Sage code to compute probabilities for multiphase distributions 6.5. Mathematica code to compute probabilities for multiphase distributions Acknowledgments References