The Branching Process: A General Conceptual Framework for Addressing Current Ecological and Evolutionary Questions

1. Introduction

Branching process is not only “the beautiful theory” in the field of probability [1], but also a key methods in evolutionary studies since its formulation [2] (Figure 1). Galton asked two questions. Suppose we start with N adult males, each with a distinct surname. In each generation, each man produces k adult sons with the probability a_k, with $\sum a_k=1$. All adults produce sons independently. What is the chance of a male whose descendants will perish after t generations? This chance is the same as the expected proportion of surnames surviving t generations. The second question is tricker. Galton wanted to have a distribution. That is, after t generations, how many surnames will have m (= 0, 1, 2, …,) bearers?

The branching process, which was used by Watson to provide a partial solution to Galton’s questions, has since found applications in ecology and evolution, through the effort of Fisher [3], Haldane [4] and many others. The basic input to the branching processes is a probability mass function specifying the reproductive potential of individuals, and the output consists of extinction probability Q, extinction time T and the variance of T [5,6]. This paper focuses on the details between the input and the output, and highlights its relevance to current research problems.

2. A Branching Process with Empirical Reproductive Success

Suppose we follow the queen and the number of her daughter queens of an invasive giant Asian hornet, Vespa mandarinia, that almost succeeded in invading North America [7] or yellow-legged Asian hornet, V. velutina that successfully invaded France and other European countries [8]. This example is also relevant to the release of animals to rescue a declining population [9] or to create a game population for hunting. Will the individuals in the new environment go extinct? How long will it take to go extinct?

Of 100 empirically monitored hornet queens, 50 did not produce any reproducing daughter queens, 30 produced one reproducing queen and 20 produced two reproducing queens. Thus, a queen produces 0, 1, and 2 daughter queens with the corresponding probabilities 0.5, 0.3 and 0.2, respectively. The probability mass function (PMF), which is often called reproduction law in branching processes, is then:

$$\begin{gathered} \mathit{f}\left(\mathit{x}=0\right)=0.5 \\ \mathit{f}\left(\mathit{x}=1\right)=0.3 \\ \mathit{f}\left(\mathit{x}=2\right)=0.2 \end{gathered}$$

(1)

We assume that this PMF does not change with time in a branching process and is applicable to every replicating queen. However, one could have branching processes to accommodate queens with different reproducing potentials.

A branching process and its analysis are anchored on a probability generating function (PGF) derived from the PFM. PGF is defined as $G_X\left(s\right)={\sum }_{k=0}^{Maxk}f\left(k\right)s^k$, where k is the value that X can take, and the subscript X can often be omitted without confusion. Given the distribution specified in Equation (1), the PGF is therefore

$$\mathit{G}\left(\mathit{s}\right)=0.5{\mathit{s}}^0+0.3{\mathit{s}}^1+0.2{\mathit{s}}^2=0.5+0.3\mathit{s}+0.2{\mathit{s}}^2$$

(2)

Note that $G_X\left(1\right)=1$. The mean and variance of X based on the PGF are

$${\mathit{E}\left(\mathit{X}\right)=\mathit{G}\mathbf{\text{'}}}_{\mathit{X}}\left(1\right)=\mathit{\mu }$$

(3)

$${\mathit{V}\mathit{a}\mathit{r}\left(\mathit{X}\right)=\mathit{G}\mathbf{\text{'}}\mathbf{\text{'}}}_{\mathit{X}}\left(1\right)+{\mathit{G}}_{\mathit{X}}^{\mathbf{\text{'}}}\left(1\right)-{\left[{\mathit{G}}_{\mathit{X}}^{\mathbf{\text{'}}}\left(1\right)\right]}^2$$

(4)

For the PGF in Equation (2), we have

$${\mathit{G}\mathbf{\text{'}}}_{\mathit{X}}\left(\mathit{s}\right)=0.3+0.4\mathit{s};{\mathit{G}}_{\mathit{X}}^{\mathbf{\text{'}}\mathbf{\text{'}}}\left(\mathit{s}\right)=0.4$$

(5)

$$\mathit{E}\left(\mathit{X}\right)={\mathit{G}}_{\mathit{X}}^{\mathbf{\text{'}}}\left(1\right)=0.3+0.4=0.7$$

(6)

$${\mathit{V}\mathit{a}\mathit{r}\left(\mathit{X}\right)=\mathit{G}\mathbf{\text{'}}\mathbf{\text{'}}}_{\mathit{X}}\left(1\right)+{\mathit{G}}_{\mathit{X}}^{\mathbf{\text{'}}}\left(1\right)-{\left[{\mathit{G}}_{\mathit{X}}^{\mathbf{\text{'}}}\left(1\right)\right]}^2=0.4+0.7-0.7^2=0.61$$

(7)

In ecological and evolutionary studies, a stable wild-type population with μ = E(X) = 1 is often implicitly assumed.

2.1. Extinction Probability Q

What is the probability that an individual who produces offspring according to Equation (1) will leave no descendants in the future? The extinction probability Q was derived [2] as the smallest non-negative solution to the following fixed-point equation:

$$\mathit{Q}=\mathit{G}\left(\mathit{Q}\right)$$

(8)

From Equation (2), we have

$$\mathit{Q}=\mathit{G}\left(\mathit{Q}\right)=0.5+0.3{\mathit{Q}}^1+0.2{\mathit{Q}}^2$$

(9)

$$0.2{\mathit{Q}}^2-0.7\mathit{Q}+0.5=0$$

(10)

This quadratic equation has two roots, $r_1=2.5,r_2=1$. Therefore, Q = 1, i.e., any population started from an individual who reproduces according to Equation (1) will certainly go extinct. Fisher [3] recognized the dependence of Q on $\mu$. If $\mu \le 1$, the mutant will ultimately die out. If $\mu >1$, then there is a non-zero chance of survival. Such findings led to the subsequent classification of subcritical ($\mu <1)$, critical ($\mu =1$) and supercritical ($\mu >1$) conditions [5,6]. Under subcritical and critical conditions, $Q=1$. Under supercritical conditions, $Q$ < 1.

The calculation above allows us to conclude that, even if the invasive hornet were lucky to produce multiple queens in certain years, it will certainly go extinct. However, how long will it take for the invasive hornet to go extinct?

2.2. Extinction Time T and Var(T)

A more difficult problem than the extinction probability is the extinction time. The individual (and its descendants) could go extinct at time t₁, t₂, t₃, etc., with their corresponding probabilities, so this extinction time should be characterized as a distribution.

Designate population size at time t as Z_t. and the extinction probability of the population by generation t as $q_t=\mathrm{Pr}(Z_t=0)$. We start the process with a single individual at time 0, so ${\mathit{Z}}_0=1$. We can compute $q_t$ by iteration using the PGF in Equation (9):

$${\mathit{q}}_{\mathit{t}+1}=\mathit{G}\left({\mathit{q}}_{\mathit{t}}\right)$$

(11)

For the first few generations:

$${\mathit{q}}_1=\mathit{G}\left({\mathit{q}}_0\right)=\mathit{G}\left(0\right)=0.5$$

(12)

$${\mathit{q}}_2=\mathit{G}\left({\mathit{q}}_1\right)=\mathit{G}\left(0.5\right)=0.5+0.3\times 0.5+0.2\times 0.5^2=0.7$$

(13)

$${\mathit{q}}_3=\mathit{G}\left({\mathit{q}}_2\right)=\mathit{G}\left(0.7\right)=0.808$$

(14)

$${\mathit{q}}_4=\mathit{G}\left({\mathit{q}}_3\right)\approx 0.876$$

(15)

The functions above introduce us to a mathematical term called composition, e.g., $q_4$ is a composition of the function G to itself 4 times. A more descriptive term is iterate [5], i.e., $q_4$ is an iterate of $q_3$ or the 4th iterate of $q_0$. As t increases, $q_t$ will increase monotonously towards Q. Under subcritical conditions, Q = 1, so $q_t$ will increase monotonously towards 1.

What is the probability that extinction happens exactly at time t, i.e., $\mathbf{Pr}\left(\mathit{T}=\mathit{t}\right)$? This $\mathbf{Pr}\left(\mathit{T}=\mathit{t}\right)$ is simply $q_t-q_{t-1}$. Thus

$$\mathbf{Pr}\left(\mathit{T}=1\right)={\mathit{q}}_1-{\mathit{q}}_0=0.5-0=0.5$$

(16)

$$\mathbf{Pr}\left(\mathit{T}=2\right)={\mathit{q}}_2-{\mathit{q}}_1=0.7-0.5=0.2$$

(17)

$$\mathbf{Pr}\left(\mathit{T}=3\right)={\mathit{q}}_3-{\mathit{q}}_2=0.808-0.7=0.108$$

(18)

$$\mathbf{Pr}\left(\mathit{T}=4\right)={\mathit{q}}_4-{\mathit{q}}_3=0.876-0.808=0.068$$

(19)

We can calculate T by the following equation:

$$\mathit{T}={\sum }_{\mathit{t}=1}^{\mathbf{\infty }}\mathit{t}\mathbf{Pr}(\mathit{T}=\mathit{t})=1\times 0.5+2\times 0.2+\dots =2.3954$$

(20)

This equation works because ${\sum }_{T=1}^{\infty }\mathrm{Pr}(T)=1$. If the extinction probability $Q<1$, then we should use the weighted average:

$$\mathit{T}={\displaystyle \frac{{\sum }_{\mathit{t}=1}^{\mathit{\infty }}\mathit{t}\mathbf{Pr}(\mathit{T}=\mathit{t})}{\mathit{Q}}}={\displaystyle \frac{{\sum }_{\mathit{t}=1}^{\mathit{\infty }}[\mathit{t}\left({\mathit{q}}_{\mathit{t}}-{\mathit{q}}_{\mathit{t}-1})\right]}{{\sum }_{\mathit{t}=1}^{\mathit{\infty }}({\mathit{q}}_{\mathit{t}}-{\mathit{q}}_{\mathit{t}-1})}}$$

(21)

An alternative way of calculating T under subcritical conditions is

$$\mathit{T}={\sum }_{\mathit{t}=1}^{\mathbf{\infty }}(1-{\mathit{q}}_{\mathit{t}})=\left(1-0\right)+\left(1-0.5\right)+\left(1-0.7\right)+\dots =2.3954$$

(22)

The variance of T can be computed by the general equation:

$$\mathit{V}\mathit{a}\mathit{r}\left(\mathit{T}\right)=\mathit{E}\left({\mathit{T}}^2\right)-{\left[\mathit{E}\left(\mathit{T}\right)\right]}^2$$

(23)

We already know $E\left(T\right)=2.3954$, so ${\left[E\left(T\right)\right]}^2=5.7379$. The first term can be calculated in either of two ways below.

$$\mathit{E}\left({\mathit{T}}^2\right)={\displaystyle \frac{{\sum }_{\mathit{t}=1}^{\mathbf{\infty }}{\mathit{t}}^2\mathbf{Pr}(\mathit{T}=\mathit{t})}{{\sum }_{\mathit{t}=1}^{\mathbf{\infty }}\mathbf{Pr}(\mathit{T}=\mathit{t})}}={\sum }_{\mathit{t}=1}^{\mathbf{\infty }}{\mathit{t}}^2\mathbf{Pr}(\mathit{T}=\mathit{t})=1^2\bullet 0.5+2^2\bullet 0.2+\dots \approx 10.7818$$

(24)

$$\mathit{E}\left({\mathit{T}}^2\right)={\sum }_{\mathit{t}=0}^{\mathbf{\infty }}\left(2\mathit{t}+1\right)(1-{\mathit{q}}_{\mathit{t}})\approx 10.7818(\mathit{u}\mathit{n}\mathit{d}\mathit{e}\mathit{r} \mathit{s}\mathit{u}\mathit{b}\mathit{c}\mathit{r}\mathit{i}\mathit{t}\mathit{i}\mathit{c}\mathit{l}\mathit{a} \mathit{c}\mathit{o}\mathit{n}\mathit{d}\mathit{i}\mathit{t}\mathit{i}\mathit{o}\mathit{n})$$

(25)

$$\mathit{V}\mathit{a}\mathit{r}\left(\mathit{T}\right)=\mathit{E}\left({\mathit{T}}^2\right)-{\left[\mathit{E}\left(\mathit{T}\right)\right]}^2=10.7818-5.7378=5.0439$$

(26)

Equation (24) is generally applicable. Equation (25) can be used only under subcritical conditions. Equation (22) and Equation (25) can be used only under subcritical conditions. In contrast, Equation (21) and Equation (24) are general and consequently can be used in subcritical, critical and supercritical conditions.

This calculation of T and its variance allows us to have more informative predictions on the fate of the invasive hornet. Not only will the invasive hornet will go extinct, but it has a probability of 0.95 of going extinct in $2.3954+1.96\times \sqrt{5.0439}$ generations. The Asian hornets, the number of generations is equivalent to the number of years.

The calculation above is summarized in Table 1, where t is the generation, $q_t$ in the second column is the probability of extinction by generation t, $q_{t+1}-q_t$ in the third column is the probability of extinction that occurs exactly in generation t, i.e., $\mathrm{Pr}(T=t)$. Note that ($q_t-q_{t-1}$) is the generation-specific extinction probability, whereas $q_t$ is the cumulative probability of extinction. T can be calculated either as the weighted average of t (weighted by $\mathrm{Pr}(T=t)$, i.e., the third column in Table 1) as shown in Equation (20), or as the summation of $1-q_t$ in the fourth column in Table 1, as shown in Equation (22). The last column is for calculating $E\left(T^2\right)$ needed for computing Var(T).

What would be the extinction time T and its Var(T) under supercritical conditions? T can be calculated the same way by using Equation (21), i.e., it is conditional on the extinction probability. Var(T) can also be calculated as before:

$$\mathit{V}\mathit{a}\mathit{r}\left(\mathit{T}\right)=\mathit{E}\left({\mathit{T}}^2\right)-{\left[\mathit{E}\left(\mathit{T}\right)\right]}^2={\displaystyle \frac{{\sum }_{\mathit{t}=1}^{\mathbf{\infty }}{\mathit{t}}^2\mathbf{Pr}\left(\mathit{T}=\mathit{t}\right)}{{\sum }_{\mathit{t}=1}^{\mathbf{\infty }}\mathbf{Pr}\left(\mathit{T}=\mathit{t}\right)}}-{\mathit{T}}^2$$

(27)

For example, in the supercritical case above with $\mathit{E}\left(\mathit{X}\right)=1.1$, we have $T\approx 5.229,Var\left(T\right)=43.826$. Note that T under the supercritical condition could be small when $\mathit{E}\left(\mathit{X}\right)$ is large. This is because such a population will either have bad luck and go extinct quickly (i.e., a small T), or increase in population size and never go extinct. For example, if the probabilities of producing 0, 1, or 2 offspring per generation are 0.1, 0.4, and 0.5, respectively, then $\mathit{E}\left(X\right)=1.5$. The resulting $T=2.159,Var\left(T\right)=3.026$.

3. A branching Process with Reproductive Success Following a Poisson Distribution

Both Fisher [3] and Haldane [4] used the Poisson distribution as the reproduction law to model reproductive success in their use of the branching process, i.e., the probability of an individual producing 0, 1, ..k, .. offspring is $\mathit{f}\left(\mathit{k}\right)={\mathit{e}}^{-\mathit{\lambda }}{\displaystyle \frac{{\mathit{\lambda }}^{\mathit{k}}}{\mathit{k}!}}$. This is also used in epidemiological studies. Let X be the number of individuals who would become infectious “descendants” of an infected individual (before the individual either dies of the disease or clears the infectious pathogen). The PGF and the associated mean and variance are

$$\mathit{G}\left(\mathit{s}\right)={\sum }_{\mathit{k}=0}^{\mathit{M}\mathit{a}\mathit{x}\mathit{k}}\mathit{f}\left(\mathit{k}\right){\mathit{s}}^{\mathit{k}}={\mathit{e}}^{-\mathit{\lambda }}+{\mathit{e}}^{-\mathit{\lambda }}\mathit{\lambda }\mathit{s}+{\mathit{e}}^{-\mathit{\lambda }}{\displaystyle \frac{{\mathit{\lambda }}^2}{2!}}{\mathit{s}}^2+{\mathit{e}}^{-\mathit{\lambda }}{\displaystyle \frac{{\mathit{\lambda }}^3}{3!}}{\mathit{s}}^3+\dots$$

(28)

$$\mathit{G}\left(\mathit{s}\right)={\mathit{e}}^{-\mathit{\lambda }}\left(1+\mathit{\lambda }\mathit{s}+{\displaystyle \frac{{\left(\mathit{\lambda }\mathit{s}\right)}^2}{2!}}+{\displaystyle \frac{{\left(\mathit{\lambda }\mathit{s}\right)}^3}{3!}}+\dots \right)={\mathit{e}}^{-\mathit{\lambda }}{\mathit{e}}^{\mathit{\lambda }\mathit{s}}={\mathit{e}}^{\mathit{\lambda }(\mathit{s}-1)}$$

(29)

$$\mathit{G}\mathbf{\text{'}}\left(\mathit{s}\right)={\mathit{e}}^{-\mathit{\lambda }}\mathit{\lambda }+2{\mathit{e}}^{-\mathit{\lambda }}{\displaystyle \frac{{\mathit{\lambda }}^2}{2!}}\mathit{s}+3{\mathit{e}}^{-\mathit{\lambda }}{\displaystyle \frac{{\mathit{\lambda }}^3}{3!}}{\mathit{s}}^2+\dots +\mathit{n}{\mathit{e}}^{-\mathit{\lambda }}{\displaystyle \frac{{\mathit{\lambda }}^{\mathit{n}}}{\mathit{n}!}}{\mathit{s}}^{\mathit{n}-1}+\dots$$

(30)

$$\mathit{G}\mathbf{\text{'}}\mathbf{\text{'}}\left(\mathit{s}\right)=2{\mathit{e}}^{-\mathit{\lambda }}{\displaystyle \frac{{\mathit{\lambda }}^2}{2!}}+6{\mathit{e}}^{-\mathit{\lambda }}{\displaystyle \frac{{\mathit{\lambda }}^3}{3!}}{\mathit{s}}^1+\dots +\mathit{n}\left(\mathit{n}-1\right){\mathit{e}}^{-\mathit{\lambda }}{\displaystyle \frac{{\mathit{\lambda }}^{\mathit{n}}}{\mathit{n}!}}{\mathit{s}}^{\mathit{n}-2}+\dots$$

(31)

$$\mathit{E}\left(\mathit{X}\right)=\mathit{G}\mathbf{\text{'}}\left(1\right)={\mathit{e}}^{-\mathit{\lambda }}\mathit{\lambda }\left(1+\mathit{\lambda }\mathit{s}+{\displaystyle \frac{{\left(\mathit{\lambda }\mathit{s}\right)}^2}{2!}}+\dots \right)=\mathit{\lambda }$$

(32)

$$\mathit{G}\mathbf{\text{'}}\mathbf{\text{'}}\left(1\right)={\mathit{e}}^{-\mathit{\lambda }}{\mathit{\lambda }}^2\left(1+\mathit{\lambda }+{\displaystyle \frac{{\mathit{\lambda }}^2}{2!}}+{\displaystyle \frac{{\mathit{\lambda }}^3}{3!}}+\dots \right)={\mathit{e}}^{-\mathit{\lambda }}{\mathit{\lambda }}^2{\mathit{e}}^{\mathit{\lambda }}={\mathit{\lambda }}^2$$

(33)

$$\mathit{V}\mathit{a}\mathit{r}\left(\mathit{X}\right)=\mathit{G}\mathbf{\text{'}}\mathbf{\text{'}}\left(1\right)+\mathit{G}\mathbf{\text{'}}\left(1\right)-{\left[\mathit{G}\mathbf{\text{'}}\left(1\right)\right]}^2={\mathit{\lambda }}^2+\mathit{\lambda }-{\mathit{\lambda }}^2=\mathit{\lambda }$$

(34)

3.1. Extinction Probability Q

Obviously, the disease propagation depends on $\mathit{\lambda }$. If $\lambda \le 1$, then disease will die out. If $\lambda >1$, then the number of descendants will increase. We will first address the same question of extinction probability Q. From what we have learned, we need to solve $Q=G\left(Q\right)$ and find the smallest non-negative root.

From Equation (29), we have

$$\mathit{Q}=\mathit{G}\left(\mathit{Q}\right)={\mathit{e}}^{\mathit{\lambda }\left(\mathit{Q}-1\right)}$$

(35)

In the critical case when $E\left(X\right)=\lambda =1$, then $\mathit{Q}={\mathit{e}}^{\left(\mathit{Q}-1\right)}$, and the only root is $Q=1$, i.e., extinction is certain. This also implies that all cases with $\lambda <1$ should also have $Q=1$ because extinction would be even more likely with smaller $\lambda$. Thus, $Q=1$ for all $\lambda \le 1$.

What would be the extinction probability under supercritical conditions with $\lambda >1$? Haldane [4] and Felsenstein [10] illustrated the solution with an $\lambda$ only slightly larger than 1 (i.e., a small selection coefficient s = $\lambda -1$) so that one can expand the right-hand side of Equation (35) as a power series and drop high-order terms. This would not work with larger $\lambda$. A better way to compute Q under supercritical conditions is to use the following equation

$$\mathit{Q}=-{\displaystyle \frac{1}{\mathit{\lambda }}}{\mathit{W}}_0\left(-\mathit{\lambda }{\mathit{e}}^{-\mathit{\lambda }}\right)$$

(36)

where W₀ is the principal branch of the Lambert W function [11]. $W_0\left(z\right)$ is defined as

$$\mathit{z}={\mathit{W}}_0\left(\mathit{z}\right){\mathit{e}}^{{\mathit{W}}_0\left(\mathit{z}\right)}$$

(37)

It can be calculated in many numerical ways [11], e.g., by using the ‘optim’ function in R or the Solver function in EXCEL. Table 2 lists Q values given different $\lambda$ values. In population genetics where $\lambda$ is the absolute fitness of a mutant, the selection coefficient s = $\lambda$ – 1 (against an implicitly assumed stable wile-type population).

Epidemiologists have documented the basic reproduction number (R₀, which is equivalent to $\lambda$ in Table 2) for various viral diseases. For COVID-19, R₀ varies from 2 to 7 [12] which corresponds to very small Q.

The calculation above is also applicable to many evolutionary scenarios, such as the fate of a mutant in a haploid population or of a parthenogenetic individual arising from a sexual population. Imagine a large, stable and genetically homogenous haploid population in which everyone is expected to produce one offspring to replace itself ($\lambda =1)$. A new mutant arises with $\lambda =1.01$, i.e., it is fitter than other individuals in the population with a fitness differential $s=0.01$. Even such a beneficial mutation would have an extinction probability $Q=0.98034$ (Table 2). Its fixation probability is $(1-Q)=1-0.98034\approx 0.02=2s$. For $s=0.001,Q=0.998002,1-Q=0.002=2s$. It is generally true that, for a small s, the fixation probability is 2s when the reproductive success is specified by a Poisson distribution.

While our illustration here is applicable to a haploid population, the conclusion is also true for diploid populations [4,10]. In a large population of wild-type individuals with AA genotypes, each individual produces 0, 1, 2, … k offspring with corresponding probabilities $p_0,p_1,p_2,\dots ,p_k$, respectively. What is Q for a new B allele that creates an AB heterozygote? If all alleles are neutral, then the probability of fixation of a new allele is simply its allele frequency, i.e., 1/(2N) for a new allele, where N is the population size of the diploid population. Thus, $Q=1-1/\left(2N\right)$. What would be the extinction probability and extinction time when the B allele is advantageous?

The two previous examples revealed that, if an advantageous individual or allele is to go extinct because of bad luck, the extinction happens quickly. If the new beneficial allele B overcomes the bad luck and becomes frequent, then it tends to achieve fixation. This allows us to assume that, during the process of the allele going extinct, its frequencies are low, and it therefore exists in heterozygotes instead of BB homozygotes. For convenience, we may also assume that the B allele is dominant or codominant so that its beneficial effect will manifest in heterozygotes. Given these conditions, the extinction probability can be written in the form that we are already familiar with:

$$\mathit{Q}={\mathit{p}}_0+{\mathit{p}}_1\mathit{Q}+{\mathit{p}}_2{\mathit{Q}}^2+{\mathit{p}}_3{\mathit{Q}}^3+\dots +{\mathit{p}}_{\mathit{k}}{\mathit{Q}}^{\mathit{k}}+\dots$$

(38)

If we further assume that $p_i$ follows the Poisson distribution, then $Q=e^{\lambda (Q-1)}$ as we have derived before in Equation (35). The solution for Q and numeric estimates for T and Var(T) have already been derived above. The assumption of independent propagation of the new allele B in the population is reflected in those terms with $Q^2,Q^3,Q^4,\dots$ In other words, all B alleles exist in AB heterozygotes and independently have the same extinction probability as the first AB heterozygote created by the mutation. With Equation (36) we can readily compute Q, T, and Var(T) for any $\lambda$.

The calculation is also applicable to the extinction probability and extinction time for a new parthenogenetic individual in a diploid sexual population. Such a parthenogenetic individual is supposed to have a two-fold advantage in fitness [13,14], i.e., with $\lambda =2$ given that an average sexual individual’s reproductive output is $\lambda =1$. The extinction probability for this new parthenogenetic individual is $Q=0.20319$ (Table 2).

3.2. Extinction Time

Let us first consider the subcritical condition with $\lambda =0.7$. This way, we can compare the results to those in the first example where $E\left(X\right)=0.7,Var\left(X\right)=0.61$. For a Poisson distribution, $E\left(X\right)=Var\left(X\right)=\lambda =0.7$. Thus, we expect the mean and variance of the extinction time to be similar to those in our first example. Again, here are the first few ${\mathit{q}}_{\mathit{t}}$ values.

$${\mathit{q}}_1=\mathit{G}\left({\mathit{q}}_0\right)=\mathit{G}\left(0\right)={\mathit{e}}^{-\mathit{\lambda }}\approx 0.4966$$

(39)

$${\mathit{q}}_2=\mathit{G}\left({\mathit{q}}_1\right)=\mathit{G}\left({\mathit{e}}^{-\mathit{\lambda }}\right)\approx 0.7030$$

(40)

$${\mathit{q}}_3=\mathit{G}\left({\mathit{q}}_2\right)\approx 0.8123$$

(41)

$${\mathit{q}}_4=\mathit{G}\left({\mathit{q}}_3\right)\approx 0.8769$$

(42)

In all subcritical cases, $q_t$ will increase monotonously towards 1. The probability that extinction happens exactly at time t, i.e., $\mathbf{Pr}\left(\mathit{T}=\mathit{t}\right)$ can also be calculated as before:

$$\mathbf{Pr}\left(\mathit{T}=1\right)={\mathit{q}}_1-{\mathit{q}}_0=0.4966$$

(43)

$$\mathbf{Pr}\left(\mathit{T}=2\right)={\mathit{q}}_2-{\mathit{q}}_1=0.2064$$

(44)

$$\mathbf{Pr}\left(\mathit{T}=3\right)={\mathit{q}}_3-{\mathit{q}}_2=0.808-0.7=0.1093$$

(45)

$$\mathbf{Pr}\left(\mathit{T}=4\right)={\mathit{q}}_4-{\mathit{q}}_3=0.876-0.808=0.0646$$

(46)

Given the distribution, we can calculate T with the two formulae that we have used before for subcritical conditions:

$$\mathit{T}={\sum }_{\mathit{t}=1}^{\mathbf{\infty }}\mathit{t}\mathbf{Pr}(\mathit{T}=\mathit{t})=1\times 0.4966+2\times 0.2064+\dots \approx 2.3762$$

(47)

$$\mathit{T}={\sum }_{\mathit{t}=0}^{\mathbf{\infty }}(1-{\mathit{q}}_{\mathit{t}})=\left(1-0\right)+\left(1-0.4966\right)+\left(1-0.7030\right)+\dots \approx 2.3762$$

(48)

What would be the extinction time T and its Var(T) under supercritical conditions? T can be calculated the same way by using Equation (21). Var(T) can also be calculated with Equation (27). For example, if $\lambda =1.1$, then $T=4.102,Var\left(T\right)=33.550$. Table 3 lists the q_t and (q_t − q_t-1) values needed to computing T and Var(T) for $\lambda =1.1$. We already have listed $Q=0.82391$ for $\lambda =1.1$ in Table 2, which is consistent with the q_t column in Table 3 that comes very close to this number when t reaches 1000 generations.

In our numeric illustration of the supercritical condition with $\lambda =1.1$, T is longer than the subcritical condition with $\lambda =0.7$. However, it is possible for a smaller T under supercritical conditions. For example, if $\lambda =2$, then $T=1.545,Var\left(T\right)=0.953$. This is because the extinction time is conditional on extinction. An individual that reproduces with $\lambda =2$ will either have bad luck and go extinct quickly (i.e., a small T), or increase in population size rapidly and never go extinct.

4. Come Back to the Galton Questions

Galton asked two questions. N adult men each have a unique surname. In each generation, each man produces k adult sons with probability a_k. Galton limited k to a maximum of five. The PGF for addressing Galton’s two questions is similar to our first example:

$$\begin{gathered} \mathit{G}\left(\mathit{s}\right)={\mathit{a}}_0+{\mathit{a}}_1\mathit{s}+{\mathit{a}}_2{\mathit{s}}^2+{\mathit{a}}_3{\mathit{s}}^3+{\mathit{a}}_4{\mathit{s}}^4+{\mathit{a}}_5{\mathit{s}}^5 \\ {\mathit{G}}^{\mathbf{\text{'}}}\left(\mathit{s}\right)={\mathit{a}}_1+2{\mathit{a}}_2\mathit{s}+3{\mathit{a}}_3{\mathit{s}}^2+4{\mathit{a}}_4{\mathit{s}}^3+5{\mathit{a}}_5{\mathit{s}}^4 \\ {\mathit{G}}^{\mathbf{\text{'}}\mathbf{\text{'}}}\left(\mathit{s}\right)=2{\mathit{a}}_2+6{\mathit{a}}_3\mathit{s}+12{\mathit{a}}_4{\mathit{s}}^2+20{\mathit{a}}_5{\mathit{s}}^4 \\ \mathit{E}\left(\mathit{X}\right)={\mathit{G}}^{\mathbf{\text{'}}}\left(1\right)={\mathit{a}}_1+2{\mathit{a}}_2+3{\mathit{a}}_3+4{\mathit{a}}_4+5{\mathit{a}}_5 \\ \mathit{V}\mathit{a}\mathit{r}\left(\mathit{X}\right)={\mathit{G}}^{\mathbf{\text{'}}\mathbf{\text{'}}}\left(1\right)+{\mathit{G}}^{\mathbf{\text{'}}}\left(1\right)-{\left[{\mathit{G}}^{\mathbf{\text{'}}}\left(1\right)\right]}^2 \end{gathered}$$

(49)

Suppose we have$a_0=0.1,a_1=0.3,a_2=0.25,a_3=0.15,a_4=0.1,a_5=0.1$. The mean number of sons produced by a man (n) and the variance of n are

$${\mathit{E}\left(\mathit{X}\right)=\mathit{G}}^{\mathbf{\text{'}}}\left(1\right)=2.15\mathit{V}\mathit{a}\mathit{r}\left(\mathit{X}\right)={\mathit{G}}^{\mathbf{\text{'}}\mathbf{\text{'}}}\left(1\right)+{\mathit{G}}^{\mathbf{\text{'}}}\left(1\right)-{\left[{\mathit{G}}^{\mathbf{\text{'}}}\left(1\right)\right]}^2=2.1275$$

(50)

This is a supercritical condition, so the population size will more than double itself each generation. The extinction probability $Q<1$. Solving the equation $Q=G\left(Q\right)$ gives us $Q=0.1519416038$. The probability of extinction by generation t is $q_t$ which we have calculated recursively in previous examples. In mathematics, this $q_t$ is calculated as the t-fold composition of G or the tth iterate of G;

$${\mathit{q}}_{\mathit{t}}={\mathit{G}}^{\left(\mathit{t}\right)}\left(0\right)$$

(51)

This function does not have a closed form, which is why we resorted to computing $q_t$ recursively in previous examples. For the first few generations:

$$\begin{gathered} {\mathit{q}}_0=0 \\ {\mathit{q}}_1=\mathit{G}\left({\mathit{q}}_0\right)=\mathit{G}\left(0\right)={\mathit{a}}_0=0.1 \\ {\mathit{q}}_2=\mathit{G}\left({\mathit{q}}_1\right)=\mathit{G}\left(0.1\right)=0.132661 \\ {\mathit{q}}_3=\mathit{G}\left({\mathit{q}}_2\right)=.1445833203 \end{gathered}$$

(52)

Thus, after one generation, 10% of the N surnames are lost. At generation 23, the proportion of surnames lost has already reached Q (=0.1519416038). After that, all remaining surnames would already have escaped bad luck and will survive forever, if the earth were infinitely large. The extinction time T is only 1.547654, with Var(T)=0.916671. This means that ~15.2% of lost surnames must have been lost in early generations, otherwise a few generations of rapid increase in population size (more than doubling in each generation) would secure their permanent existence.

Thus, to address Galton’s first question, as long as $a_i$ values and t are given, it is easy to recursively compute $q_t$, which is the proportion of surnames that have gone extinct by generation t (Table 4). For example, about 14.9% of the surnames are expected to perish by generation 4 under the extremely favorable condition with $\mu =2.15$.

Galton’s second question is trickier. At generation t, how many of those starting N surnames are carried by 0, 1, 2, …, men? Our illustrative example has $a_0=0.1,a_1=0.3,a_2=0.25,a_3=0.15,a_4=0.1,a_5=0.1$. Suppose we have N = 10000 men each bearing a unique surname at time t = 0. At time t = 1, we expect 1000 ($=N\bullet a_0$) surnames to perish. We expect 3000, 2500, 1500, 1000 and 1000 surnames, respectively, to have 1, 2, 3, 4 and 5 bearers, respectively (Figure 2A). At time t = 2, the proportion of perished surnames is expected to increase from 10% to 13.277% (Figure 2B). The distribution is also wider (Figure 2B), because there are some lucky surnames that might have as many as 25 ($=5\times 5$) bearers after the second generation. For example, the man with a surname Clark might produce 5 sons in the first generation, and these sons each produce 5 sons in the second generation. Now the surname Clark has 25 bearers. However, the probability of having such good luck is small, i.e., the product of $a_5$ from the first generation multiplied by $a_5^5$ in the second generation (i.e., each of the five sons independently produces five sons in the second generation. There are about 1159 surnames with two bearers after the second generation (Figure 2B). After the third generation, the largest possible number of bearers of a surname is 125 (=5³) (Figure 2C). However, the probability for such an event to occur is 10⁻³¹. The total probability of having a surname with 50 or more bearers is only 3.24⁻⁴. Note that this is under extremely favorable conditions with $\mu =2.15$. In a stable population with $\mu =1$ (e.g., when $a_0=0.45,a_1=0.35,a_2=0.15,a_3=0.06,a_4=0.03,a_5=0.01$), then the surnames lost would be 45% after the first generation, 64% after the second generation, and 78% after the third generation. The total probability for surnames with 10 or more bearers is only 2.97⁻³ after three generations.

5. Recent Extensions and Applications

The largest number of illustrations and potential biological application of the branching process is available in a book [15]. Recent studies have extended the classical branching process in several ways. The first is parameter estimation. Equation (1) listed three parameter values, 0.5, 0.3 and 0.2 that defines the reproduction law needed for the branching process. The λ parameter in the Poisson distribution is also a parameter that can be estimated from empirical data. Such estimation can be done either nonparametrically or by maximum likelihood or Bayesian approaches [16], the last being particularly relevant for small populations with limited empirical data.

Second, the branching process itself can be generalized. In our illustration, the same reproduction law is applied to all individuals. However, parental reproductive success could be correlated with offspring reproductive success. Such a scenario implies that some lineages are less likely to go extinct than others. For example, successful mothers tend to have successful daughters and so on. One can modify the reproduction law to accommodate this parent-offspring correlation in reproductive success [17]. Also, my illustrative examples are for populations with discrete generations. Under certain circumstances, the branching process can also accommodate populations with overlapping generations [17].

The theory of multi-type branching processes has been developed a long time ago [5,6], but only recently been applied to solve practical biological problems on the consequences of dormancy on population demography. A haploid population was implicitly assumed. The empirical challenge is to determine what environmental factors trigger the switch between the active and dormant state and whether the probability of switching between the two states depends quantitatively on the degree of environmental harshness.