Branching Processes with Mating, Part I

Up to now we have ignored that in species with sexual reproduction changes in population size depend on the formation of couples. Sometimes this may be justified. For instance, if mating is polygynous (one male may have many females), and the number of males is sufficiently large, the total number of females in a generation determines future population size. Another example is the number of initially rare mutants in an evolutionary process, for which, provided the residence population is large enough, pair formation is not a limiting factor. In many instances, however, mating is an important factor that cannot be dismissed. Bisexual Galton–Watson processes take mating into account explicitly. Unfortunately, this complicates life considerably, and few mathematical results are available.

So-called bisexual branching processes – a somewhat unfortunate name for processes with mating (a.k.a. two sex branching processes) – were introduced by Daley (1968). In them it is convenient to represent population size by the number of couples rather than the number of individuals. Denote the number of couples in the \(n\)th generation by \(Z_n\). The process starts with an initial number of couples, \(Z_0\), just before the first reproduction period. Each couple has random numbers of female and male offspring. Thus, the offspring of a couple is represented by a pair of numbers, which are independently distributed for different couples, with an identical distribution that is constant over generations. In the simplest case it is assumed that, conditional on the total number of offspring, the offspring sex is determined at random (like flipping a, possibly biased, coin). Another model of sex determination assumes that the numbers of male and female offspring are independent with a possibly different distribution. The two mechanisms are equivalent for Poisson-distributed numbers of male and female offspring, since the sum of two independent Poisson variables is Poisson again. As in the Galton–Watson process, generations are non-overlapping.

The most realistic scenario is to form pairs according to a stochastic process. In existing theory, however, pair formation is assumed to be deterministic. Thus, the formation of couples is determined by a mating function \(\zeta (x, y)\), which specifies the total number of couples formed if there are \(x\) females and \(y\) males. If mating is strictly monogamous and if we assume that the population is well-mixed, so that all possible pairs are formed, \(\zeta (x, y)\) equals the minimum of \(x\) and \(y\). In a polygamous mating system, in which each male may inseminate up to \(d\) females, \(\zeta (x, y)\) is equal to the minimum of \(x\) and \(dy\). Another, less common example, is unilateral promiscuous mating, where \(\zeta (x, y)\) equals \(x\) times \(\min (1, y)\). All these examples belong to the class of so-called common sense mating functions, which satisfy three natural conditions, i.e.,

Furthermore, these functions are super-additive (Hull 1982), which means that for all \(x_1\), \(x_2\), \(y_1\), \(y_2=0,1,2,\ldots\) \[ \zeta (x_1+x_2, y_1+y_2) \ge \zeta (x_1, y_1)+\zeta (x_2, y_2). \] Note Super-additivity of the mating function implies that, if the population is divided into two subgroups that each form couples, the total number of couples is at most equal to that of the undivided population.

We number the couples in the nth generation in an arbitrary way and let \(X_k\) and \(Y_k\) denote the numbers of female and male offspring, respectively, of the \(k\)th couple. Conditionally, on the number of couples in the \((n−1)\)th generation, \(Z_{n−1}\), the number of couples in the next generation equals \[ Z_n=\zeta \left(\sum_{k=1}^{Z_{n-1}} X_k, \sum_{k=1}^{Z_{n-1}} Y_k\right). \] A process is extinct when \(Z_n\) equals zero. Bruss (1984) pointed out that the growth of a population is determined by the average unit reproduction means \[ \mu_j=\frac{1}{j}E(Z_n\ |\ Z_{n-1}=j), \qquad j\ge 1. \]
Note In the standard Galton-Watson process the expectation in the equation above equals \(j\) times \(\mu\), so the average unit reproduction means \(\mu_j\) equal the offspring mean \(\mu\).

Daley et al. (1986) showed that for common sense, superadditive mating functions, as \(j\) tends to infinity, \(\mu_j\) converges to a limit \(\mu_\infty\) with \(\mu_j\le \mu_\infty\) for all \(j\), and that extinction is certain (irrespective of the initial population size) if this number is less than or equal to one. If \(mu_\infty>1\), there is a positive survival probability provided the initial population size is large enough. If, however, the mating function and reproduction functions are such that, for all \(i\) \[ P(Z_n>i\ |\ Z_{n-1}=i)>0, \] then \(\mu_\infty>1\) implies a ositive ssurvival probability for any initial positive ppopulation size.

From the above we can conclude \[ E(Z_n\ | Z_{n-1})=\mu_{Z_{n-1}} Z_{n-1}\le \mu_\infty Z_{n-1}, \] and then \[ E(Z_n)\le \mu_\infty E(Z_{n-1}\le \ldots \le \mu_\infty^n E(Z_{n-1}). \] Thus, the expected population size does not necessarily change exponentially, as in the asexual branching process model already considered. Instead, decline may be more rapid if \(\mu_\infty \le 1\), and growth may be slower in the case of \(\mu_\infty >1\). The growth rate depends on characteristics of the mating function and offspring distributions.

Example Consider a small population of swans. Since these birds are strictly monogamous, we can use a model with mating function \(\zeta (x, y)=\min(x, y)\). We assume that each couple has a chance \(1−q\) of producing zero offspring and a chance \(q\) of producing three surviving offspring. Given that a couple has offspring, the number of females is distributed \(Bin(3,p)\), with \(p=0.5\). Denote by \(P(x,y\ |\ j)\) the conditional probability that, given the adult population size equals \(j\), the total number of female offspring equals \(x\) and the total number of males equals \(y\). Using the formula \(E(X)=\sum_{k\ge 1}P(X\ge k)\), which holds for any non-negative integer-valued \(X\), we obtain \[ E(\min(X,Y)\ |\ j)=\sum_{k=1}^{3j} P(\min (X,Y)\ge k\ |\ j)=\sum_{k=1}^{3j}\sum_{x=k}^{3j}\sum_{y=k}^{3j}P(x,y\ |\ j). \]

The probabilities \(P(x,y\ |\ j)\) can best be derived using probability generating functions. Alternatively, one can use Mathematica to do these rather tedious calculations. Figure 1 shows \(\mu_j\) as a function of \(j\) for two different values of \(q\). When \(q\) equals \(0.7\), \(\mu_\infty\) is about 0.85, so the population certainly becomes extinct. The expected number of offspring per couple is 2.1 in this case, so the expected number of offspring per individual equals 1.05 and an asexual branching process with this value of \(\mu\) would have a positive survival probability. When \(q\) equals 0.9, \(\mu_\infty\) equals approximately 1.1 and, accordingly, there is a positive chance that the population will persist. The survival probability of the population also depends on the proportion of females that is produced. When the chance of female offspring is not 1/2, the value of \(\mu_\infty\) drops, since the number of pairs that can be formed becomes lower. As shown in Figure 1, when \(q = 0.9\) the process becomes nearly critical when the chance of having daughters, \(p\), is 0.4 rather than 0.5. This illustrates that in sexually reproducing species, not only are the expected numbers of individuals important, but also their sex ratio is a critical factor in determining extinction probability. Obviously, this also depends on the mating system. When there is (slight) polygyny, for instance, we can expect that female-biased sex ratios are less risky than male-biased ones. The magnitude of such effects may be studied with bisexual branching process models.

Figure 1 Values of \(\mu_j\) for the model given in the Example, with different parameter values.