5.2. MultiType Galton-Watson processes: Generation-based framework

Consider a Galton-Watson process with two types, in which we think of offspring as occuring in discrete generations. In principle, each type may have some number of offspring of either type.

We use

\[ \xi_1(x,y) = \sum_{j=0}^\infty \sum_{k=0}^\infty p_{j,k}x^jy^k \]

to denote the offspring distribution of the first type and

\[ \xi_2(x,y) = \sum_{j=0}^\infty \sum_{k=0}^\infty q_{j,k}x^jy^k \]

to denote the offspring distribution of the second type.

5.2.1. Joint PGF

Let \((X_g,Y_g)\) represent \(X_g\) individuals of the first type and \(Y_g\) individuals of the second type at generation \(g\).

We assume that the process begins with a single individual of one of the two types. The PGF of the distribution at generation \(g\) given a single initial individual of the first type is written

\[ \Phi_g(x,y|(1,0)) = \sum_{j,k} p_{j,k}(g) x^jy^k \]

where

\[ p_{j,k}(g)=\mathbb{P}[X_g=j, Y_g=k | X_0=1, Y_0=0] \]

Similarly we write the PGF given a single individual of the second type

\[ \Phi_g(x,y|(0,1)) = \sum_{j,k} q_{j,k}(g) x^jy^k \]

where

\[ q_{j,k}(g)=\mathbb{P}[X_g=j, Y_g=k | X_0=0, Y_0=1] \]

This notation \(p\) and \(q\) overlaps with the coefficients listed for \(\xi_1\) and \(\xi_2\). We will see that in terms of the coefficients of the PGFs \(\xi_1(x,y)\) and \(\xi_2(x,y)\), we have \(p_{j,k}(1) = p_{j,k}\) and \(q_{j,k}(1)=q_{j,k}\).

We use \(\vec{\Phi}_g(x,y)\) to denote \(\left(\Phi_g(x,y|(1,0)), \quad\Phi_g(x,y|(0,1))\right)\). We find that

\[\begin{align*} \Phi_1(x,y|(1,0)) &= \xi_1(x,y)\\ \Phi_2(x,y|(0,1)) &= \xi_2(x,y) \end{align*}\]

So

\[ \vec{\Phi}_1(x,y) = (\xi_1(x,y), \xi_2(x,y)) \]

As commented above, this gives that \(p_{j,k}(1) = p_{j,k}\) and \(q_{j,k}(1) = q_{j,k}\).

To go to higher generations we have two options.

5.2.1.1. Backward approach

Our first method resembles the Backward Kolmogorov Equation derivation. We consider first \(\Phi_{g+1}(x,y|(1,0))\). We know that the PGF for generation \(1\) from this initial state is \(\Phi_{1}(x,y|(1,0)) = \xi_1(x,y)\). That is, the joint distribution of \(X_1, Y_1\), the numbers of each type at generation \(1\) has PGF \(\Phi_1(x,y|(1,0))\). Each individual of type \(1\) produces a number of descendents \(g\) generations later (i.e., at generation g+1) whose PGF is \(\Phi_{g}(x,y|(1,0))\). Similarly each individual of type \(2\) produces a number whose PGF is \(\Phi_g(x,y|(0,1))\). Taking the randomly-stopped sum of offspring from \(X_1\) individuals of type \(1\) and \(Y_1\) individuals of type \(2\) yields the PGF

\[\begin{align*} \Phi_{g+1}(x,y|(1,0)) &= \xi_1\left(\Phi_{g}(x,y|(1,0)),\quad \Phi_{g}(x,y|(0,1))\right)\\ &= \xi_1(\vec{\Phi}_{g-1}(x,y)) \end{align*}\]

and similarly \(\Phi_g(x,y|(0,1)) = \xi_2(\vec{\Phi}_{g-1}(x,y))\). Thus we can recursively define

\[ \vec{\Phi}_g(x,y) = \left(\xi_1(\vec{\Phi}_{g-1}(x,y)), \quad \xi_2(\vec{\Phi}_{g-1}(x,y))\right) \]

with \(\vec{\Phi}_0(x,y) = (x,y)\).

In particular, this means that

Theorem 5.3 (PGF of size distribution of two-type Galton-Watson process at generation \(g\) (Backward version))

Given a two-type Galton-Watson process, with offspring distributions \(\xi_1(x,y)\) and \(\xi_2(x,y)\), the PGFs for the distribution at generation \(g\) can be found by recursively solving

\[ \vec{\Phi}_{g+1}(x,y) = \left(\xi_1(\vec{\Phi}_{g}(x,y)),\quad \xi_2(\vec{\Phi}_{g}(x,y))\right) \]

with

\[ \vec{\Phi}_0(x,y) = (x,y) \]

Example 5.1 (Joint distribution of number infected or recovered at generation \(g\))

Assume that an outbreak begins with a single infected individual in a large enough population that we can assume the Galton-Watson assumptions hold at least through generation \(g\). We let \(p_{j,k}(g)\) denote the probability of \(j\) active infections and \(k\) recoveries by generation \(g\).

Assume also that the number of infections for the next generation caused by an infected individual has PGF \(\mu(x)\). The infected individual recovers at the next generation, producing one recovered individual. So the PGF for the offspring distribution counting infections and recoveries at the next generations is \(\xi_1(x,y) = y\mu(x)\).

A recovered individual at one generation remains a recovered individual at the next generation, so \(\xi_2(x,y)=y\).

Thus \(\vec{\Phi}_g(x,y) = \left(\xi_1(\vec{\Phi}_{g-1}(x,y)),\quad \mu_2(\vec{\Phi}_{g-1}(x,y))\right)\). Looking at the first few generations, taking \(\xi_1(x,y) = y\mu(x)\) and \(\xi_2(x,y)=y\) we have

\[\begin{align*} \vec{\Phi}_0(x,y) &= (x,y)\\ \vec{\Phi}_1(x,y) &= \left(\xi_1(\vec{\Phi}_0(x,y)),\quad \xi_2(\vec{\Phi}_0(x,y))\right)\\ &= (\xi_1(x,y),\xi_2(x,y))\\ &= (y\mu(x),y)\\ \vec{\Phi}_2(x,y) &= \left(\xi_1(\vec{\Phi}_1(x,y)),\quad \xi_2(\vec{\Phi}_1(x,y))\right)\\ &= \left(\xi_1(y\mu(x),y),\quad \xi_2(y \mu(x),y)\right)\\ &= (y \mu(y\mu(x)), \quad y)\\ & \: \vdots \end{align*}\]

So \(\Phi_g(x,y|(1,0)) = y \mu\left(\Phi_{g-1}(x,y|(1,0))\right)\).

If we want the PGF of the number recovered at generation \(g\), we have \(\Phi_g(1,y|(1,0)) = y \mu(\Phi_{g-1}(1,y|(1,0)))= y\mu(y\mu(\cdots y\mu(1) \cdots))\).

5.2.1.2. Forward Approach

The second method resembles the Forward Kolmogorov Equation derivation. We consider first \(\Phi_{g+1}(x,y|(1,0))\). We know that the PGF for generation \(g\) is \(\Phi_g(x,y|(1,0))\). This means that the joint distribution of \(X_g,Y_g\), the numbers of each type at generation \(g\) has PGF \(\Phi_g(x,y|(1,0))\). Each of these produces a random number at generation \(g+1\) whose PGFs are \(\xi_1(x,y)\) and \(\xi_2(x,y)\) for type \(1\) or type \(2\) respectively. Thus the PGF for the randomly-stopped sum of offspring from \(X_g\) individuals of type \(1\) and \(Y_g\) of type \(2\) is

\[ \Phi_{g+1}(x,y|(1,0)) = \Phi_g(\xi_1(x,y), \xi_2(x,y) | (1,0)) \]

Symmetrically appling this to \(\Phi_g(x,y|(0,1))\) we conclude

Theorem 5.4 (PGF of size distribution of two-type Galton-Watson process at generation \(g\) (Forward version))

Given a two-type Galton-Watson process, with offspring distributions \(\xi_1(x,y)\) and \(\xi_2(x,y)\), the PGFs for the distribution at generation \(g\) can be found by recursively solving

\[ \vec{\Phi}_{g+1}(x,y) = \vec{\Phi}_g(\xi_1(x,y), \xi_2(x,y)) \]

with

\[ \vec{\Phi}_0(x,y) = (x,y) \]

Example 5.2 (Joint distribution of number infected or recovered at generation \(g\))

Assume that an outbreak begins with a single infected individual in a large enough population that we can assume the Galton-Watson assumptions hold at least through generation \(g\). We let \(p_{j,k}(g)\) denote the probability of \(j\) active infections and \(k\) recoveries by generation \(g\).

Assume also that the number of infections for the next generation caused by an infected individual has PGF \(\mu(x)\). The infected individual recovers at the next generation, producing one recovered individual. So the PGF for the offspring distribution counting infections and recoveries at the next generations is \(\xi_1(x,y) = y\mu(x)\).

A recovered individual at one generation remains a recovered individual at the next generation, so \(\xi_2(x,y)=y\).

Thus \(\vec{\Phi}_g(x,y) = \vec{\Phi}_{g-1}(\xi_1(x,y),\xi_2(x,y)) = \vec{\Phi}_{g-1}(y\mu(x),y)\) with \(\Phi_0(x,y) = (x,y)\). We conclude

\[\begin{align*} \vec{\Phi}_0(x,y) &= (x,y)\\ \vec{\Phi}_1(x,y) &= \vec{\Phi}_0\left(\xi_1(x,y),\quad\xi_2(x,y)\right)\\ &= \vec{\Phi}_0(y\mu(x), y)\\ &= (y \mu(x), y)\\ \vec{\Phi}_2(x,y) &= \vec{\Phi}_1\left(\xi_1(x,y),\quad\xi_2(x,y)\right)\\ &= \vec{\Phi}_1(y \mu(x),y)\\ &= (y \mu(y\mu(x)), y)\\ & \:\vdots \end{align*}\]

So \(\Phi_{g}(x,y|(1,0)) = \Phi_{g-1}(\mu(x),y|(1,0)) = y\mu(y\mu(\cdots y\mu(x,y) \cdots))\).

5.2.2. Extinction Probability

In this context, extinction can have multiple meanings. We may mean extinction of either type or extinction of both types.

We begin by thinking of extinction of both types. We assume, without loss of generatlity that the outbreak begins with a single individual of type \(1\). If \(\Phi_g(x,y|(1,0)) = \sum_{j,k} p_{j,k} x^jy^k\) then the probability that both types are extinct by generation \(g\) is

\[\begin{align*} \mathbb{P}[\text{extinction of both types by generation $g$}|(1,0)] &= p_{0,0}(g)\\ &= \sum_{j,k} p_{j,k}0^j0^k\\ &= \Phi_g(0,0|(1,0)) \end{align*}\]

and similarly given the initial individual being of type \(2\).

Using the Backward framework, we have

\[ \vec{\Phi}_g(0,0)= \left(\xi_1(\vec{\Phi}_{g-1}(0,0)), \quad \xi_2(\vec{\Phi}_{g-1}(0,0))\right) \]

Setting \(\vec{\alpha}(g) = \vec{\Phi}_g(0,0)\), this can be solved iteratively:

\[ \vec{\alpha}(g) = \left(\xi_1(\vec{\alpha}(g-1)), \quad\xi_2(\vec{\alpha}(g-1))\right) \]

with \(\vec{\alpha}(0) = (0,0)\).

Using the Forward framework, we have

\[\begin{align*} \vec{\Phi}_g(0,0) &= \vec{\Phi}_{g-1}\left(\xi_1(0,0), \xi_2(0,0)\right)\\ &= \vec{\Phi}_{g-1}(p_{0,0}, q_{0,0}) \end{align*}\]

This cannot be solved iteratively.

Theorem 5.5 (Probability of complete extinction of two-type Galton-Watson Process)

The probabilities of extinction of both types by generation \(g\), depending on initial condition is \(\vec{\alpha}(g) = (\alpha(g|(1,0)), \alpha(g|(0,1)))\). This is found by iteratively solving

\[ \vec{\alpha}(g)= \left(\xi_1(\vec{\alpha}(g-1)), \quad\xi_2(\vec{\alpha}(g-1))\right) \]

with \(\vec{\alpha}(0) = (0,0)\).

The long-time extinction probability from each state solve \(\vec{\alpha} = (\xi_1(\vec{\alpha}), \vec{\xi}_2(\vec{\alpha}))\).

Note - I want to come back to this following bit about 1 type going extinct — but I want to wait until I see if it’s really a useful result for later applications. We can really just restict our attention to the reduced PGF setting the other variable to be 1

If we just want the extinction of one type, we restrict our attention to the case where the other type has only itself for offspring. That is, we do not consider extinction of just one type if it could reappear as offspring of the other type. So we assume that if we are interested in extinction of type \(1\), we must have \(\xi_2(x,y) = \sum_{j,k} q_{j,k} x^j y^k\) is independent of \(x\). So \(q_{j,k} = 0\) for any \(j>0\). We write \(\xi_2(x,y) =\xi_2(y)= \sum_{k} q_{0,k} y^k\)

In this case we want to find \(\vec{\Phi}_g(0,y)\). Again the backward version is more useful. We know that if the initial individual is of type \(2\) then type \(1\) is immediately extinct and remains extinct. Thus setting \(\alpha(g;1)\) to be the extinction probability of type \(1\) if the initial individual is of type \(1\) we have

\[\begin{align*} \alpha &= \Phi_g(0,y|(1,0))\\ &= \xi_1(\alpha(g-1;1),1) \end{align*}\]

with \(\alpha(0;1)=0\).

Theorem 5.6 (Probability of extinction of one type in a two-type Galton-Watson Process)

If \(\xi_2(x,y)\) can be written \(\xi_2(y)\), then the probabilities of extinction of type \(1\) by generation \(g\) assuming an initial individual of type \(1\) is \(\alpha(g;1)\) where

\[ \alpha(g;1) = \xi_1(\alpha(g-1;1),1) \]

with \(\alpha(0;1) = 0\).

The symmetric result holds for \(\alpha(g;2)\), the probability that the second type goes extinct assuming that \(\xi_1(x,y) = \xi_1(x)\) and the initial individual is of type \(2\).

5.2.3. Expected size

We look at the matrix of expected values \(\mathbb{E}[X_g]\) and \(\mathbb{E}[Y_g]\) assuming \((X_0,Y_0)\) is \((1,0)\) and \((0,1)\). We could work this out without going through generating functions. In hindsight, that is probably easier than what I did below.

\[\begin{align*} \begin{pmatrix} \mathbb{E}[X_g|(1,0)],\quad \mathbb{E}[X_g|(0,1)]\\ \mathbb{E}[Y_g|(1,0)],\quad \mathbb{E}[Y_g|(0,1)] \end{pmatrix} =& \begin{pmatrix} \left[\frac{\partial}{\partial x} \Phi_{g}(x,1 |(1,0))\right|_{x=1}& \left[\frac{\partial}{\partial x} \Phi_{g}(x,1 |(0,1))\right|_{x=1}\\ \left[\frac{\partial}{\partial y} \Phi_{g}(1,y|(1,0))\right|_{y=1} & \left[\frac{\partial}{\partial y} \Phi_{g}(1,y|(0,1))\right|_{y=1}\end{pmatrix}\\ % %=& %\begin{pmatrix} \left[\frac{\partial}{\partial x} \vec{\Phi}_{g-1} \left( \xi_1(x,1), \xi_2(x,1)\right)\right|_{x=1}\\ %\left[\frac{\partial}{\partial y} \vec{\Phi}_{g-1} \left( \xi_1(1,y), \xi_2(1,y)\right)\right|_{y=1}\end{pmatrix}\\ \end{align*}\]

We now focus just on the first column, giving the probabilities assuming the first individual is of type \(1\). We note that the second column follows through symmetry. We use the Forward version of the equations for which \(\Phi_g(x,y|(1,0)) = \Phi_{g-1}(\xi_1(x,y), \xi_2(x,y)|(1,0))\)

\[\begin{align*} \begin{pmatrix}\mathbb{E}[X_g|(1,0)] \\ \mathbb{E}[Y_g|(1,0)]\end{pmatrix} % =& % \begin{pmatrix} \left[\frac{\partial}{\partial x} \Phi_{g-1}(\xi_1(x,1),\xi_2(x,1)|(1,0))\right|_{x=1}\\ \left[\frac{\partial}{\partial y} \Phi_{g-1}(\xi_1(1,y),\xi_2(1,y)|(1,0))\right|_{y=1} \end{pmatrix}\\ % =& \begin{pmatrix} \left[\left(\frac{\partial}{\partial \xi_1},\frac{\partial}{\partial \xi_2} \right) \Phi_{g-1} \left( \xi_1(x,1), \xi_2(x,1)|(1,0) \right) \cdot \frac{\partial}{\partial x} \left(\xi_1(x,1), \xi_2(x,1)\right)\right|_{x=1}\\ \left[\left(\frac{\partial}{\partial \xi_1},\frac{\partial}{\partial \xi_2} \right) \Phi_{g-1} \left( \xi_1(1,y), \xi_2(1,y)|(1,0) \right) \cdot \frac{\partial}{\partial y} \left(\xi_1(1,y), \xi_2(1,y)\right)\right|_{y=1}\end{pmatrix}\\ % =& \begin{pmatrix} \left[\left(\frac{\partial}{\partial \xi_1},\frac{\partial}{\partial \xi_2} \right) \Phi_{g-1} \left( \xi_1, \xi_2|(1,0) \right)\right|_{\xi_1,\xi_2=1,1} \cdot \left[\frac{\partial}{\partial x} \left(\xi_1(x,1), \xi_2(x,1)\right)\right|_{x=1}\\ \left[\left(\frac{\partial}{\partial \xi_1},\frac{\partial}{\partial \xi_2} \right) \Phi_{g-1} \left( \xi_1, \xi_2|(1,0) \right)\right|_{\xi_1,\xi_2=1,1} \cdot \left[\frac{\partial}{\partial y} \left(\xi_1(1,y), \xi_2(1,y)\right)\right|_{y=1}\end{pmatrix}\\ % =&\begin{pmatrix}\left[\frac{\partial}{\partial x} \xi_1(x,1)\right|_{x=1} & \left[\frac{\partial}{\partial x} \xi_2(x,1)\right|_{x=1}\\ \left[\frac{\partial}{\partial y} \xi_1(1,y)\right|_{y=1}& \left[\frac{\partial}{\partial y}\xi_2(1,y)\right|_{y=1}\end{pmatrix} \begin{pmatrix}\left[\frac{\partial}{\partial \xi_1} \Phi_{g-1}(\xi_1,1|(1,0))\right|_{\xi_1=1}\\ \left[\frac{\partial}{\partial \xi_2} \Phi_{g-1}(1,\xi_2|(1,0))\right|_{\xi_2=1}\end{pmatrix} \\ % =&\begin{pmatrix}\mathbb{E}[X|(1,0)]& \mathbb{E}[X|(0,1)]\\ \mathbb{E}[Y|(1,0)]& \mathbb{E}[Y|(0,1)]\end{pmatrix} \begin{pmatrix}\mathbb{E}[X_{g-1}|(1,0)]\\ \mathbb{E}[Y_{g-1}|(1,0)]\end{pmatrix}\\ =& \begin{pmatrix}\mathbb{E}[X|(1,0)]& \mathbb{E}[X|(0,1)]\\ \mathbb{E}[Y|(1,0)]& \mathbb{E}[Y|(0,1)]\end{pmatrix}^g \begin{pmatrix}\mathbb{E}[X_{0}|(1,0)]\\ \mathbb{E}[Y_{0}|(1,0)]\end{pmatrix}\\ =& \begin{pmatrix}\mathbb{E}[X|(1,0)]& \mathbb{E}[X|(0,1)]\\ \mathbb{E}[Y|(1,0)]& \mathbb{E}[Y|(0,1)]\end{pmatrix}^g \begin{pmatrix}1\\0\end{pmatrix} \end{align*}\]

Where \(\mathbb{E}[X|(1,0)]\) is the expected number of direct \(X\) offspring of one \(X\) individual, etc.

So including the second column, we have

\[\begin{split} \begin{pmatrix} \mathbb{E}[X_g|(1,0)],\quad \mathbb{E}[X_g|(0,1)]\\ \mathbb{E}[Y_g|(1,0)],\quad \mathbb{E}[Y_g|(0,1)] \end{pmatrix}= \begin{pmatrix}\mathbb{E}[X|(1,0)]& \mathbb{E}[X|(0,1)]\\ \mathbb{E}[Y|(1,0)]& \mathbb{E}[Y|(0,1)]\end{pmatrix}^g \begin{pmatrix}1 & 0\\0 &1\end{pmatrix} \end{split}\]