Stochastic Models of Stem Cells and Their Descendants under Different
Criticality Assumptions

Nam H Nguyen111CONTACT<EMAIL_ADDRESS>and Marek Kimmel

Department of Statistics, Rice University, Houston, TX, USA

###### Abstract

We study time continuous branching processes with exponentially distributed
lifetimes, with two types of cells that proliferate according to binary
fission. A range of possible system dynamics are considered, each of which is
characterized by the mutation rate of the original cells and the survival
probability of the altered cells’ progeny. For each system, we derive a
closed-form expression for the joint probability generating function of cell
counts, and perform asymptotic analysis on the behaviors of the cell
population with particular focus on probability of extinction. Part of our
results confirms known properties of branching processes using a different
approach while other are original. While the model is best suited for modeling
the fate of differentiating stem cells, we discuss other scenarios in which
these system dynamics may be applicable in real life. We also discuss the
history of the subject.

Keywords: stochastic processes, population dynamics, multi-type branching
processes, exact representation

##  1 Introduction

All cells existing in muticellular organisms originate from a single
fertilized egg. In embryonic development, the descendants of this single cell
migrate and differentiate to create diverse tissues and organs. Even in a
mature organism, the functional cells which perform a range of tasks have to
renew, since they are “used up”, and their building materials such as nucleic
acids, proteins and others are recycled or disposed of. In most cases, this
self-renewal system has a hierarchical structure, with the self-renewing stem
cells ($A$-cells) at the top, which divide into two progeny, each of which may
remain a stem cell or change (differentiate) into a commited cell ($B$-cell).
It is known that stem cells are highly protected and that stem cells usually
have exactly 2 progeny. In addition to this, under physiological conditions,
the number of stem $A$-cells remains roughly constant, until the organism
begins to age. The commitment is irreversible, and commited $B$-cells
proliferate so that a constant flux of their descendants is ensured. The count
of the committed $B$-cells remains roughly constant, too. In the presence of a
constant flux of differentiated $A$-cells, this is only possible if the
fraction of $B$-cells proceeding to self-renewal after division is less than
$1$. If it is greater than or equal to $1$, then the system increases in
volume (cell count), starting from the $B$-cells. Therefore, the basic system
that is interesting under physiological conditions is, in the branching
process language, a critical process of ${A}$-cells feeding into the sub-
critical process of ${B}$-cells. A sufficient condition for this to happen is
that each $A$-cell progeny becomes (on average) a $B$-cell with probability
equal to $1/2$, and each $B$-cell progeny remains (on average) a $B$-cell with
probability less than $1/2$.

However, in the early development phase, during aging, or under other
conditions in which differentiation of $B$-cells is altered, the sub-
criticality of $B$-cells may not be ensured. Less likely, but not
inconceivably, the $A$-cells might be sub- or super-critical (albeit just
slightly). Hence, we will consider a range of possibilities. The process can
be described as follows

$\displaystyle A\rightarrow AA\qquad$ $\displaystyle\text{with probability}\
(1-\alpha)^{2}$ $\displaystyle A\rightarrow AB\qquad$ $\displaystyle\text{with
probability}\ 2\alpha(1-\alpha)$ $\displaystyle A\rightarrow BB\qquad$
$\displaystyle\text{with probability}\ \alpha^{2}$ $\displaystyle B\rightarrow
BB\qquad$ $\displaystyle\text{with probability}\ p^{2}$ (1.0.1) $\displaystyle
B\rightarrow B\qquad$ $\displaystyle\text{with probability}\ 2pq$
$\displaystyle B\rightarrow\emptyset\qquad$ $\displaystyle\text{with
probability}\ q^{2}$

where $\alpha\in[0,1]$, $p\in[0,1]$, and $q=1-p$. In addition, we assume that
the lifetimes of $A$-cells and $B$-cells are exponentially distributed with
parameters $\lambda_{A}$ and $\lambda_{B}$ respectively.

The paper is organized as follows. In Section 2, we introduce notations that
will be used for the rest of the paper. In Section 3, we mention useful
properties of special functions, which will be used in the mathematical
derivations. In Section 4.1, we perform exploratory analysis to understand how
cell counts behave on average over time. Section 4.2 contains the main text of
the paper, where we derive closed-form expressions for the joint probability
generating functions that fully describe the population dynamics at any given
time. Based on the explicit results, we study the asymptotic extinction
probabilities of the system under different combinations of criticalities in
Section 4.3. In Section 5, we analyze the results and provide literature
examples of diseases which might be, perhaps metaphorically, due to all types
of deviations from the critical $A$-cells, sub-critical $B$-cells stereotype.
Section 6 is the conclusion of the paper.

##  2 Notations

Let $F_{A}(x,y,t)$ be the joint probability generating function (PGF) for the
number of cells at time $t$ given that the process is initiated by a single
cell of type $A$. That is, let $Z_{A}(t)$ and $Z_{B}(t)$ be the number of
$A$-cells and $B$-cells at time $t$ respectively. Then, $F_{A}(x,y,t)$ is
defined as follows

$\displaystyle F_{A}(x,y,t)=E\big{[}x^{Z_{A}(t)}y^{Z_{B}(t)}\big{]}$

Since we start with one single cell of type $A$, the initial condition for
$F_{A}$ is $F_{A}(x,y,0)=x$. Similarly, $F_{B}(x,y,t)$ is the joint PGF given
that the process is initiated by a single cell of type $B$. Note that cells of
type $A$ mutate irreversibly to cells of type $B$. Hence, $F_{B}(x,y,t)$ is
independent of $x$, and hence can be written more succintly as $F_{B}(y,t)$.
For consistency of notation, we will write $F_{B}(x,y,t)$ for the rest of the
paper. For $F_{B}$, we start with one cell of type $B$, so the initial
condition is $F_{B}(x,y,0)=y$.

Let $h_{A}(x,y)$ and $h_{B}(x,y)$ be the progeny PGFs for cells of type $A$
and cells of type $B$ respectively. From (1), it follows that

$h_{A}(x,y)=[(1-\alpha)x+\alpha y]^{2}$ (2.0.1) $h_{B}(x,y)=(p+qy)^{2}$
(2.0.2)

##  3 Special functions

Throughout the paper, we will make use of various special functions. Let
$I(a,z)$ be the modified Bessel function of order $a$ evaluated at $z$. The
first derivative of $I(a,z)$ is given by

$\frac{d}{dz}I(a,z)=\frac{1}{2}(I(a-1,z)+I(a+1,z))$ (3.0.1)

The following recursive formulae is also useful to simplify the results

$I(a-1,z)-I(a+1,z)=\frac{2a}{z}I(a,z)$ (3.0.2)

For asymptotic analysis, we need the following large-$z$ approximation of
$I(a,z)$

$I(a,z)\sim\frac{1}{\sqrt{2\pi z}}e^{z}$ (3.0.3)

Let $M(a,b,z)$ and $W(a,b,z)$ be the Whittaker functions (first kind and
second kind, respectively) with parameters $a$ and $b$ evaluated at $z$. The
first derivatives are given by

$\displaystyle
M^{\prime}(a,b,z)=\Big{(}\frac{1}{2}-\frac{a}{z}\Big{)}M(a,b,z)+\frac{\frac{1}{2}+a+b}{z}M(a+1,b,z)$
$\displaystyle
W^{\prime}(a,b,z)=\Big{(}\frac{1}{2}-\frac{a}{z}\Big{)}W(a,b,z)-\frac{1}{z}W(a+1,b,z)$
(3.0.4)

Let ${\vphantom{F}}_{1}\kern-0.5ptF_{1}(a,b,z)$ and $U(a,b,z)$ denote the
confluent hypergeometric functions (first kind and second kind, respectively)
with parameter $a$ and $b$ evaluated at $z$. There is a simple algebraic
relationship between the Whittaker functions and the confluent hypergeometric
functions

$\displaystyle
M(a,b,z)=e^{-z/2}z^{a+1/2}{\vphantom{F}}_{1}\kern-0.5ptF_{1}\bigg{(}b-a+\frac{1}{2},1+2b,z\bigg{)}$
$\displaystyle
W(a,b,z)=e^{-z/2}z^{b+1/2}U\bigg{(}b-a+\frac{1}{2},1+2b,z\bigg{)}$ (3.0.5)

The large-$z$ behavior of ${\vphantom{F}}_{1}\kern-0.5ptF_{1}(a,b,z)$ and
$U(a,b,z)$ is given by

$\displaystyle{\vphantom{F}}_{1}\kern-0.5ptF_{1}(a,b,z)=\frac{\Gamma(b)}{\Gamma(a)}e^{z}z^{a-b}(1+O(|z|^{-1}))$
$\displaystyle U(a,b,z)=z^{-a}(1+O(|z|^{-1}))$ (3.0.6)

Lastly, let ${\vphantom{F}}_{2}\kern-0.5ptF_{1}(a,b,c,z)$ be the Gaussian
hypergeometric function with parameters $a$, $b$ and $c$ evaluated at $z$. The
function is defined by the following series expansion

${\vphantom{F}}_{2}\kern-0.5ptF_{1}(a,b,c,z)=\frac{\Gamma(c)}{\Gamma(a)\Gamma(b)}\sum_{n=0}^{\infty}\frac{\Gamma(a+n)\Gamma(b+n)}{\Gamma(c+n)}\frac{z^{n}}{n!}$
(3.0.7)

with convergence guaranteed within the unit circle $|z|=1$. The first
derivative of ${\vphantom{F}}_{2}\kern-0.5ptF_{1}(a,b,c,z)$ is given by

${\vphantom{F}}_{2}\kern-0.5ptF_{1}^{\prime}(a,b,c,z)=\frac{ab}{c}\;{\vphantom{F}}_{2}\kern-0.5ptF_{1}(1+a,1+b,1+c,z)$
(3.0.8)

Note that the parameters of the special functions above can take complex
values. Details on the properties of these special functions can be found in
[1, 12]

##  4 Results

###  4.1 Analysis of expectations

Let $E_{A}(t)$ be the expected number of cells of type $A$ at time $t$.
Similarly, let $E_{B}(t)$ be the expected number of cells of type $B$ at time
$t$. From (1), we have that $E_{A}$ satisfies

$\frac{dE_{A}}{dt}=\lambda_{A}(1-2\alpha)E_{A}$

Under the initial condition $E_{A}(0)=1$, the solution is
$E_{A}(t)=e^{\lambda_{A}(1-2\alpha)t}$. Note that $E_{A}(t)$ is identically
$1$ when $\alpha=1/2$, which is expected from a critical process. From the
dynamics of cells of type $B$ in (1), we have that $E_{B}(t)$ must satisfy the
following first-order linear differential equation

$\frac{dE_{B}}{dt}=2\lambda_{A}(1-\alpha)e^{\lambda_{A}(1-2\alpha)t}+\lambda_{B}(q^{2}-p^{2})E_{B}$

At time $t=0$, we have no cells of type $B$. Therefore, the initial condition
is simply $E_{B}(0)=0$. If all cells are critical (i.e., $\alpha=1/2$ and
$p=q=1/2$; we call this case “bi-critical”), the equation above simplifies,
and the general solution is $E_{B}(t)=\lambda_{A}t$. If the system is not bi-
critical, we obtain the general solution

$E_{B}(t)=\frac{2\lambda_{A}(1-\alpha)}{\lambda_{A}(1-2\alpha)-\lambda_{B}(q^{2}-p^{2})}e^{\lambda_{A}(1-2\alpha)t}+Ce^{\lambda_{B}(q^{2}-p^{2})t}$

Under the initial condition $E_{B}(0)=0$, the particular solution is given by

$E_{B}(t)=\frac{2\lambda_{A}(1-\alpha)}{\lambda_{A}(1-2\alpha)-\lambda_{B}(q^{2}-p^{2})}\Big{(}e^{\lambda_{A}(1-2\alpha)t}-e^{\lambda_{B}(q^{2}-p^{2})t}\Big{)}$

###  4.2 Explicit probability generating function solution

To derive the explicit solution, we consider the following backward equations
[17]:

$\frac{dF_{A}}{dt}=-\lambda_{A}F_{A}+\lambda_{A}h_{A}(F_{A},F_{B})$ (4.2.1)
$\frac{dF_{B}}{dt}=-\lambda_{B}F_{B}+\lambda_{B}h_{B}(F_{A},F_{B})$ (4.2.2)

To solve this system of backward equations, we will solve Equation (4.2.2) for
$F_{B}$ first and then substitute it into Equation (4.2.1) to obtain $F_{A}$.
The approach that we employ is largely inspired by [2, 3], which also
concerned closed-form expressions for the PGFs of two-type time continuous
branching processes with exponentially distributed lifetimes, but with
different dynamics from ours. A two-stage model of carcinogenesis, which
involved three cell types, was addressed in [7] by reducing the partial
differential equation satisfied by the PGF to the hypergeometric differential
equation of Gauss.

######  Theorem 1.

Given cells of type $A$ and $B$ that proliferate according to dynamics (1),
the joint PGF $F_{A}(x,y,t)$ under different criticalities is given as
follows.

(a) (Bi-critical) Let $\mu=\frac{4\lambda_{A}}{\lambda_{B}}$.

$\displaystyle F_{A}(w)=\frac{-4}{\mu}\bigg{[}$
$\displaystyle\frac{\mu}{2\sqrt{w}}\frac{\big{(}\frac{-\theta}{\mu\sqrt{w}}I_{0}^{-}(w)+I_{1}^{-}(w)\big{)}+C\big{(}\frac{\theta}{\mu\sqrt{w}}I_{0}^{+}(w)+I_{1}^{+}(w)\big{)}}{I_{0}^{-}(w)+CI_{0}^{+}(w)}$
$\displaystyle+\frac{1}{2w}-\frac{\mu(1+w)}{4w}\bigg{]}$ (4.2.3)

where $w=\frac{\lambda_{B}t(1-y)/4+1}{1-y}$, $\theta=\sqrt{1-\mu}$ is a
constant defined in terms of the model parameters, and $C=C(x)$ is a constant
that is determined by the initial condition. Furthermore,

$\displaystyle I_{0}^{-}(w)=I(-\theta,\mu\sqrt{w});\
I_{0}^{+}(w)=I(\theta,\mu\sqrt{w});$ $\displaystyle
I_{1}^{-}(w)=I(-\theta+1,\mu\sqrt{w});\ I_{1}^{+}(w)=I(\theta+1,\mu\sqrt{w})$
(4.2.4)

(b) (Non-critical A-cells and critical B-cells) Let
$\mu=\frac{4\lambda_{A}}{\lambda_{B}}$.

$\displaystyle F_{A}(w)=\frac{-1}{\mu(1-\alpha)^{2}}\bigg{[}$
$\displaystyle\frac{\theta}{2}-\frac{\theta_{1}}{w}+\frac{1}{w}\frac{\big{(}\frac{1}{2}+\theta_{1}+\theta_{2}\big{)}M_{1}(w)-CW_{1}(w)}{M_{0}(w)+CW_{0}(w)}$
$\displaystyle-\frac{\mu(w+2\alpha(1-\alpha)(1-w))}{2w}\bigg{]}$ (4.2.5)

where $w=\frac{\lambda_{B}t(1-y)/4+1}{1-y}$, and $\theta_{1}$, $\theta_{2}$
and $\theta$ are constants defined in terms of the model parameters

$\theta_{1}=-\frac{\mu(1-\alpha)\alpha}{|1-2\alpha|};\
\theta_{2}=\frac{1}{2}\sqrt{1-4\mu\alpha+4\mu\alpha^{2}};\
\theta=|1-2\alpha|\mu$

$C$ is a constant that is determined by the initial condition, and

$\displaystyle M_{0}(w)=M(\theta_{1},\theta_{2},\theta w);\quad
M_{1}(w)=M(\theta_{1}+1,\theta_{2},\theta w);$ $\displaystyle
W_{0}(w)=W(\theta_{1},\theta_{2},\theta w);\quad
W_{1}(w)=W(\theta_{1}+1,\theta_{2},\theta w)$ (4.2.6)

(c) (Non-critical B-cells) Let $\mu=\frac{\lambda_{A}}{\lambda_{B}(p-q)}$.

$\displaystyle F_{A}(z)=\frac{-z}{\mu(1-\alpha)^{2}}\bigg{\\{}$
$\displaystyle\bigg{[}\frac{1+\theta_{1}}{2}z^{-\frac{1-\theta_{1}}{2}}F_{0}^{+}(z)+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}z^{\frac{1+\theta_{1}}{2}}F_{1}^{+}(z)$
$\displaystyle+C\bigg{(}\frac{1-\theta_{1}}{2}z^{-\frac{1+\theta_{1}}{2}}F_{0}^{-}(z)+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}z^{\frac{1-\theta_{1}}{2}}F_{1}^{-}(z)\bigg{)}\bigg{]}$
$\displaystyle\bigg{/}\bigg{[}z^{\frac{1+\theta_{1}}{2}}F_{0}^{+}(z)+Cz^{\frac{1-\theta_{1}}{2}}F_{0}^{-}(z)\bigg{]}+\frac{1}{2}\frac{1+\theta_{2}}{z-q^{2}}$
$\displaystyle-\frac{1}{2}\bigg{[}\frac{1}{z}-\frac{\mu}{z}\bigg{(}2\alpha(1-\alpha)\frac{z-p^{2}}{z-q^{2}}-1\bigg{)}\bigg{]}\bigg{\\}}$
(4.2.7)

where $z=\frac{p^{2}-yq^{2}}{1-y}e^{\lambda_{B}(p-q)t}$, and $\theta_{1}$,
$\theta_{2}$ and $\theta_{3}$ are constants in terms of the model parameters

$\theta_{1}=\frac{\mu\sqrt{q^{2}-4\alpha(1-\alpha)p^{2}}}{q};\
\theta_{2}=\frac{\sqrt{q^{2}-4\alpha(1-\alpha)\mu(p-q)}}{q};\
\theta_{3}=\mu|1-2\alpha|$

$C$ is a constant that is determined by the initial condition, and

$\displaystyle
F_{0}^{+}(z)={\vphantom{F}}_{2}\kern-0.5ptF_{1}\bigg{(}\frac{1}{2}(1+\theta_{1}+\theta_{2}+\theta_{3}),\frac{1}{2}(1+\theta_{1}+\theta_{2}-\theta_{3}),1+\theta_{1},\frac{z}{q^{2}}\bigg{)};$
$\displaystyle
F_{0}^{-}(z)={\vphantom{F}}_{2}\kern-0.5ptF_{1}\bigg{(}\frac{1}{2}(1-\theta_{1}+\theta_{2}+\theta_{3}),\frac{1}{2}(1-\theta_{1}+\theta_{2}-\theta_{3}),1-\theta_{1},\frac{z}{q^{2}}\bigg{)};$
$\displaystyle
F_{1}^{+}(z)={\vphantom{F}}_{2}\kern-0.5ptF_{1}\bigg{(}1+\frac{1}{2}(1+\theta_{1}+\theta_{2}+\theta_{3}),1+\frac{1}{2}(1+\theta_{1}+\theta_{2}-\theta_{3}),2+\theta_{1},\frac{z}{q^{2}}\bigg{)};$
$\displaystyle
F_{1}^{-}(z)={\vphantom{F}}_{2}\kern-0.5ptF_{1}\bigg{(}1+\frac{1}{2}(1-\theta_{1}+\theta_{2}+\theta_{3}),1+\frac{1}{2}(1-\theta_{1}+\theta_{2}-\theta_{3}),2-\theta_{1},\frac{z}{q^{2}}\bigg{)}$
(4.2.8)

###### Proof.

(a) Under bi-criticality, we obtain the following backward equation for
$B$-cells from (2.0.2) and (4.2.2)

$\frac{dF_{B}}{dt}=\frac{\lambda_{B}}{4}(F_{B}-1)^{2}$

Noting that the initial condition is $F_{B}(x,y,0)=y$, we have

$F_{B}(x,y,t)=\frac{\lambda_{B}t(1-y)/2+2y}{\lambda_{B}t(1-y)/2+2}$ (4.2.9)

Combining (2.0.1), (4.2.1) and (4.2.9) gives us

$\frac{dF_{A}}{dt}=-\lambda_{A}F_{A}+\lambda_{A}\bigg{(}(1-\alpha)F_{A}+\alpha\frac{\lambda_{B}t(1-y)/2+2y}{\lambda_{B}t(1-y)/2+2}\bigg{)}^{2}$

We then employ the transformation
$f=F_{B}=\frac{\lambda_{B}t(1-y)/2+2y}{\lambda_{B}t(1-y)/2+2}$, which leads to

$\frac{dF_{A}}{df}=\mu\frac{(1-\alpha)^{2}}{(1-f)^{2}}F_{A}^{2}+\mu\frac{2\alpha(1-\alpha)f-1}{(1-f)^{2}}F_{A}+\mu\alpha^{2}\frac{f^{2}}{(1-f)^{2}}$

where $\mu=\frac{4\lambda_{A}}{\lambda_{B}}$ is a constant. We employ another
transformation, $w=\frac{1}{1-f}$, to simplify the equation. Note that
$\frac{dw}{df}=\frac{1}{(1-f)^{2}}=w^{2}$. Therefore, the equation becomes

$\frac{dF_{A}}{dw}=\mu(1-\alpha)^{2}F_{A}^{2}+\mu\bigg{(}2\alpha(1-\alpha)\frac{w-1}{w}-1\bigg{)}F_{A}+\mu\alpha^{2}\bigg{(}\frac{w-1}{w}\bigg{)}^{2}$

We have that $w(t=0)=\frac{1}{1-y}$. Hence, the initial condition becomes
$F_{A}\Big{(}w=\frac{1}{1-y}\Big{)}=x$. This is a non-linear Riccati
differential equation, which is not possible to solve in general. However, in
this case we can convert it to a second-order linear ODE (the Sturm-Liouville
equation) and try to solve it instead [5]. To do this, we introduce the
following short-hand notations

$A=\mu(1-\alpha)^{2};\ B=\mu\bigg{(}2\alpha(1-\alpha)\frac{w-1}{w}-1\bigg{)};\
C=\mu\alpha^{2}\bigg{(}\frac{w-1}{w}\bigg{)}^{2}$

Following [3], we perform the transformation
$F_{A}=\frac{-u^{\prime}}{Au}=-\frac{1}{A}(\log u)^{\prime}$. The equation
becomes

$u^{\prime\prime}+\gamma u^{\prime}+\beta u=0$

where the coefficients $\gamma$ and $\beta$ are given by

$\gamma=-\bigg{(}\frac{A^{\prime}}{A}+B\bigg{)}=\mu\bigg{(}1+2\alpha(1-\alpha)\frac{1-w}{w}\bigg{)};\
\beta=AC=\mu^{2}\alpha^{2}(1-\alpha)^{2}\bigg{(}\frac{1-w}{w}\bigg{)}^{2}$

To remove the first derivative, we use another quasi-linear transformation
$u=\Phi T$ such that $\Phi^{\prime}=-\gamma\Phi/2$ [3]. As we shall see, we do
not need to compute the function $\Phi$ explicitly. Finally, we arrive at an
equation whose solution can be obtained explicitly

$T^{\prime\prime}-\bigg{(}\frac{\mu^{2}(w+4\alpha(1-\alpha)(1-w))}{4w}-\frac{\mu\alpha(1-\alpha)}{w^{2}}\bigg{)}T=0$
(4.2.10)

Since we are assuming that cells of type $A$ behave critically, the equation
above reduces to

$T^{\prime\prime}-\frac{\mu}{4w}\Big{(}\mu-\frac{1}{w}\Big{)}T=0$

The closed-form solution for the ODE above can be written as a linear
combination of modified-Bessel functions of the first kind. The solution, up
to a multiplicative constant, is given by

$T(w)=\sqrt{w}[I(-\sqrt{1-\mu},\mu\sqrt{w})+CI(\sqrt{1-\mu},\mu\sqrt{w})]$

where $I(.)$ is the modified Bessel function defined in Section 2. Using the
short-hand notations (1), we can write
$T(w)=\sqrt{w}[I_{0}^{-}(w)+CI_{0}^{+}(w)]$. To obtain an explicit expression
for $F_{A}$, we reverse the calculations as follows

$F_{A}(w)=\frac{-1}{A}\bigg{(}\frac{T^{\prime}(w)}{T(w)}-\frac{\gamma}{2}\bigg{)}$
(4.2.11)

As we can see, we are interested in the ratio between $T^{\prime}(w)$ and
$T(w)$ and, therefore, the omitted multiplicative constant does not matter.
Using formulas (3.0.1) and (3.0.2), and short-hand notations (1), we arrive at
the following expression for $T^{\prime}(w)$

$T^{\prime}(w)=\frac{\mu}{2}\bigg{[}\bigg{(}\frac{-\theta}{\mu\sqrt{w}}I_{0}^{-}(w)+I_{1}^{-}(w)\bigg{)}+C\bigg{(}\frac{\theta}{\mu\sqrt{w}}I_{0}^{+}(w)+I_{1}^{+}(w)\bigg{)}\bigg{]}+\frac{1}{2w}T(w)$

where $\theta=\sqrt{1-\mu}$. From (4.2.11), given that $A=\frac{\mu}{4}$ and
$\gamma=\frac{\mu(1+w)}{2w}$ in the critical case, we find that $F_{A}(w)$ is
given by (1), where $w=\frac{1}{1-f}=\frac{\lambda_{B}t(1-y)/4+1}{1-y}$ is a
function of $t$. The initial condition $F_{A}\big{(}w=\frac{1}{1-y}\big{)}=x$
leads to the following expression for the constant $C$.

$C=\frac{\kappa
I_{0}^{-}\big{(}\frac{1}{1-y}\big{)}-\big{(}\frac{-\theta\sqrt{1-y}}{\mu}I_{0}^{-}\big{(}\frac{1}{1-y}\big{)}+I_{1}^{-}\big{(}\frac{1}{1-y}\big{)}\big{)}}{\big{(}\frac{\theta\sqrt{1-y}}{\mu}I_{0}^{+}\big{(}\frac{1}{1-y}\big{)}+I_{1}^{+}\big{(}\frac{1}{1-y}\big{)}\big{)}-\kappa
I_{0}^{+}\big{(}\frac{1}{1-y}\big{)}}$ (4.2.12)

where

$\kappa=\frac{2-y-x}{2\sqrt{1-y}}-\frac{\sqrt{1-y}}{\mu}$ (4.2.13)

∎

Proofs of the remaining parts follow the same approach, and are given in
details in the Supplement.

###  4.3 Asymptotic analysis

In this section, we will study the asymptotic properties of the systems we
consider, with particular focus on the probability of extinction of cells in
large time.

######  Theorem 2.

Suppose that cells of type $A$ and $B$ proliferate according to dynamics (1).

(a) (Bi-critical) Eventual extinction occurs with probability $1$.
Furthermore, the rate at which the population approaches extinction is
proportional to the square root of time.

(b) (Non-critical A cells and critical B cells) Eventual extinction happens
with probability $\frac{\alpha^{2}}{(1-\alpha)^{2}}$ in the super-critical
case, and with probability $1$ in the sub-critical case. In both cases, the
rate at which the probability of extinction approaches the limit is
proportional to time.

(c) (Super-critical B cells) Assuming $\theta_{1}\neq-1$, where $\theta_{1}$
is defined in part (c) of Theorem 1, eventual extinction happens with
probability
$\frac{1}{2(1-\alpha)^{2}}\Big{[}1-\frac{2\alpha(1-\alpha)p^{2}}{q^{2}}-\frac{\sqrt{q^{2}-4\alpha(1-\alpha)p^{2}}}{q}\Big{]}$.
The probability of extinction approaches this limit exponentially with respect
to time.

###### Proof.

(a) The probability of extinction of the entire population by time $t$,
$E(t)$, is given by $F_{A}(0,0,t)$. From (4.2.13), we have that, at $x=0$ and
$y=0$,

$\kappa_{0}\coloneqq\kappa|_{x=0,y=0}=1-\frac{1}{\mu}$

Therefore, from (4.2.12),

$C_{0}\coloneqq
C|_{x=0,y=0}=\frac{\kappa_{0}I_{0}^{-}(1)-\big{(}\frac{-\theta}{\mu}I_{0}^{-}(1)+I_{1}^{-}(1)\big{)}}{\big{(}\frac{\theta}{\mu}I_{0}^{+}(1)+I_{1}^{+}(1)\big{)}-\kappa_{0}I_{0}^{+}(1)}$

is a finite constant. Given that $w=\frac{\lambda_{B}t(1-y)/4+1}{1-y}$, we
have that $w_{0}\coloneqq w|_{x=0,y=0}=\frac{\lambda_{B}}{4}t+1$. Hence,
$w_{0}\to\infty$ as $t\to\infty$. Using the large-argument approximation
(3.0.3) for modified Bessel functions and the fact that $C_{0}$ is finite, we
deduce that

$\frac{\big{(}\frac{-\theta}{\mu\sqrt{w_{0}}}I_{0}^{-}(w_{0})+I_{1}^{-}(w_{0})\big{)}+C_{0}\big{(}\frac{\theta}{\mu\sqrt{w_{0}}}I_{0}^{+}(w_{0})+I_{1}^{+}(w_{0})\big{)}}{I_{0}^{-}(w_{0})+C_{0}I_{0}^{+}(w_{0})}\sim
1-\frac{\theta}{\mu\sqrt{w_{0}}}\frac{1-C_{0}}{1+C_{0}}$

for large time. Therefore, it follows that

$\displaystyle E(t)$
$\displaystyle=\frac{-4}{\mu}\bigg{[}\frac{\mu}{2\sqrt{w_{0}}}\frac{\big{(}\frac{-\theta}{\mu\sqrt{w_{0}}}I_{0}^{-}(w_{0})+I_{1}^{-}(w_{0})\big{)}+C_{0}\big{(}\frac{\theta}{\mu\sqrt{w_{0}}}I_{0}^{+}(w_{0})+I_{1}^{+}(w_{0})\big{)}}{I_{0}^{-}(w_{0})+CI_{0}^{+}(w_{0})}$
$\displaystyle\qquad\quad+\frac{1}{2w_{0}}-\frac{\mu(1+w_{0})}{4w_{0}}\bigg{]}$
$\displaystyle\sim\frac{-4}{\mu}\bigg{[}\frac{\mu}{2\sqrt{w_{0}}}\bigg{(}1-\frac{\theta}{\mu\sqrt{w_{0}}}\frac{1-C_{0}}{1+C_{0}}\bigg{)}+\frac{1}{2w_{0}}-\frac{\mu(1+w_{0})}{4w_{0}}\bigg{]}$
$\displaystyle\sim\frac{-4}{\mu}\bigg{(}\frac{\mu}{2\sqrt{w_{0}}}-\frac{\mu}{4}\bigg{)}$
$\displaystyle\sim 1-\frac{4/\sqrt{\lambda_{B}}}{\sqrt{t}}$

The theorem then follows. ∎

Certain extinction of a two-type system with double criticalities has been
proven in [4]. Hence, in the bi-critical case, our theorem agrees with
established theoretical results in literature using a different approach.
Detailed proofs of the remaining parts are given in the Supplement.

##  5 Discussion

###  5.1 Stochastic models of hematopoiesis over time

The history of stochastic models of hematopoiesis goes back to the work of
Till, McCullogh, and Siminovitch in 1963 [27], who modeled proliferation of
colony-forming units (CFU) and described the distribution of sizes of
resulting colonies using gamma distribution and simulated (using an early IBM
machine) birth-death processes. Macken and Perelson [19], in a Springer
Lecture Notes volume, proposed a series of more comprehensive models of
hematopietic system in the form of a multi-type Galton-Watson process. This
analysis was insightful, including probabilities of non-extinction
(“completion of growth”), growth rates and the first two moments of cell
counts. The results were mostly obtained using now classical theorems from
Harris’s book [14]. Many original papers on stochastic models appeared since,
including Dingli et al. [8] regarding stochastic clonal expansion of
hematopoietic stem cells, Kimmel and Corey [18] and Wojdyla et al. [31]
regarding evolution of neutropenia-related leukemias, and many others. Some of
them were reviewed in [30, 16]. A paper that features a three-type branching
process model not very unlike ours was published by Denes and Krewski [7].

Diverse aspects of the dynamics of the hematopoietic systems can be understood
using deterministic or quasi-stochastic approaches. The former involves
expected-value dynamics modeled using ordinary differential equation (ODE)
models, while the latter involves expected frequencies modeled using integral
equations or partial differential equations of transport type. A comprehensive
review of both types of approaches was carried out by Pujo-Menjouet [22]. A
recent paper by Dinh et al. [9] shows how stochasticity is important in
prediction of outcomes of treatments for leukemias (for the deterministic
model related to the latter paper, see [26]).

###  5.2 Review of biologically relevant cases

We will limit the review to the haemopoietic (blood cell production) system in
the Human and in laboratory animals such as the Mouse. The $A$-cells here are
the haemopoietic stem cells (HST), while the $B$-cells are the committed
multipotent progenitors (CMP), which have the ability to feed into more
specialized compartments, with probability $p$ for each progeny cell following
division. CMP’s have self-renewal capacity, with probability $q$ for each
progeny cell following division, but they mainly serve the role of
transitional cells [30]. Analogies with other hierarchical cell production
systems, such as the neurogenic progenitors in the brain [28], urothelial
cells in the urinary bladder [29], and others, are immediate. However, it
should be noted that, except for the neurogenesis, no other cell production
system seems to reach the sophistication level of higher mammals
haematopoiesis. Please note that we do not consider transients caused by
regulatory feedbacks, such as the interferon $\gamma$ feedback in
haematopoiesis [15], which appears when the system is under stress.

Normal adult haematopoiesis: In this case, stabilization of the mean HSC and
CMP counts is essential [23]. This is only achieved if the $A$-cells are
critical ($\alpha=1/2$) and the $B$-cells are sub-critical ($p>q$).

Normal fetal haematopoiesis: In this case, growth of the mean HSC and CMP
counts is essential [31]. This type of haematopoiesis occurs in the Human over
the last 3 months of fetal development. This may be achieved if the $A$-cells
are super-critical ($\alpha>1/2$) and the $B$-cells are sub-critical, critical
or super-critical. In each of the 3 cases, the relative proportion of HSC to
CMP will be changing differently with time. The probability of extinction in
any finite time depends on the type of proliferation of $B$-cells. Inroads
into understanding this process have been obtained using mosaic mice
methodology [10].

Normal aging haematopoiesis: In this case, the mean HSC count is decreasing
but the CMP count should be prevented from falling too fast [23]. This type of
haematopoiesis occurs in the Human at age greater than about 70. This may be
achieved if the $A$-cells are sub-critical ($\alpha>1/2$) and the $B$-cells
are super-critical ($p<q$).

Repopulation following bone marrow transplant: This case is not very
dissimilar from the fetal haematopoiesis, as the main purpose is growth of the
mean HSC and CMP counts. An interesting twist is that the transplant has to
include a mix of some HSC and much more CMP, since transplanting HSC alone
will not allow building up of haematopoiesis in a short time; growth of the
mean HSC and CMP counts is essential [25]. Our model allows ballpark
computations of this mix, given the minimum count of CMPs to be present after
a fixed time from the transplant to prevent subject’s death.

Cyclic neutropenias: Mathematically interesting cases coincide frequently with
biological reality. For example, in the bi-critical case, we expect large
oscillations, observed in diseases such as cyclic neutropenias [21]. In the
historic paper by Mackey and Glass [20], cyclic hematopoiesis is attributed to
defects in feedback governing stem cell activation, but the strong stochastic
component may be better explained by criticality.

###  5.3 Distributions of cell lifetimes

Methods, which we used to obtain our theorems, are typical for a linear birth-
and-death process [13]. As it is known, the latter is mathematically identical
to a binary fission branching process with exponential lifetimes, which
traditionally is considered not appropriate for modeling cell populations.
Indeed, in laboratory cell populations, the sojourn times in different cell
cycle phases are usually of finite duration, and in some bacteria, cell
divisions are almost synchronous, as attested by a voluminous literature of
the subject (e.g., Kimmel and Axelrod’s monograph [17]). However, the
situation is different in cancer and physiological control systems. For
example, in lung cancer, time to doubling of cell populations is of the order
of several weeks [11]. In haematopoietic stem cells, it may be even longer
[6]. Much of this time is waiting for “permission” to divide (in healthy
cells) or finishing of multiple rounds of DNA repair (in cancer cells). This
part of the cell cycle is likely to be distributed exponentially, although
measurements in vivo are difficult. In any case, models of physiological stem
cell systems based on exponentially distributed cell lifetimes (hence expected
values are described by ordinary differential equations) seem to work with
considerable precision. This was demonstrated among others by measurements and
mathematical models based on shortening of telomeres in haematopoietic stem
cells [24].

##  6 Conclusions

Two- and multi-type branching processes are minimalist models for normal cell
production systems and cancer progression in which altered clones play a major
role. We presented closed-form solutions for biologically meaningful two-type
time continuous branching processes with exponentially distributed lifetimes.
We computed the probability generating functions explicitly for all possible
combinations of cell-type-specific criticalities. The most practically
interesting case, which arises in the dynamics of stem cells, is when the
self-renewing stem cells are critical and the commited cells are sub-critical.
However, as discussed previously, other system dynamics are also relevant in
certain situations. While systems with non-critical dynamics are complex due
to the presence of the confluent hypergeometric and Whittaker functions,
solution for the bi-critical case can be expressed in terms of the relatively
simple modified Bessel functions with well-defined asymptotic behaviors. Our
asymptotic analyses not only reveal the large-time probability of extinction,
but also describe the rate at which this probability approaches the limit.
Given the closed-form solutions, other interesting asymptotic properties can
be studied by considering certain scaling limits, and require further research
in the future. It is also more realistic to model cancer progression using a
multi-type branching process. Given the already complex nature of the explicit
solutions presented in our paper, doing the same for a multi-type system may
not be feasible (however, see [7]) and require a different approach.

## Acknowledgments

Nam Nguyen acknowledges funding from CPRIT grant RP200383 (Dr. Wenyi Wang,
PI). Marek Kimmel acknowledges funding from NIH R01HL136333 and R01HL134880
grants (Dr. Katherine King, PI).

## References

  * [1] Milton Abramowitz and Irene A. Stegun.  Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables.  Dover, New York, ninth edition, 1964.
  * [2] Tibor Antal and Pavel Krapivsky.  Exact solution of a two-type branching process: Clone size distribution in cell division kinetics.  J Stat Mech, 2009.
  * [3] Tibor Antal and Pavel Krapivsky.  Exact solution of a two-type branching process: Models of tumor progression.  J Statist Mech Theory Exp, 8, 2011.
  * [4] K.B. Athreya, P.E. Ney, and P.E. Ney.  Branching Processes.  Dover Books on Mathematics. Dover Publications, 2004.
  * [5] Carl Bender and Steven Orszag.  Advanced Mathematical Methods for Scientists and Engineers: Asymptotic Methods and Perturbation Theory, volume 1.  Springer, New York, NY, 1999.
  * [6] Sandra N Catlin, Lambert Busque, Rosemary E Gale, Peter Guttorp, and Janis L Abkowitz.  The replication rate of human hematopoietic stem cells in vivo.  Blood, The Journal of the American Society of Hematology, 117(17):4460–4466, 2011.
  * [7] Josef Denes and Daniel Krewski.  An exact representation for the generating function for the moolgavkar-venzon-knudson two-stage model of carcinogenesis with stochastic stem cell growth.  Mathematical Biosciences, 131(2):185–204, 1996.
  * [8] David Dingli, Arne Traulsen, and Jorge M. Pacheco.  Stochastic dynamics of hematopoietic tumor stem cells.  Cell Cycle, 6:461–466, 2007.
  * [9] Khanh N Dinh, Roman Jaksik, Seth J Corey, and Marek Kimmel.  Predicting time to relapse in acute myeloid leukemia through stochastic modeling of minimal residual disease based on clonality data.  Computational and Systems Oncology, 1(3):e1026, 2021.
  * [10] Miguel Ganuza, Trent Hall, David Finkelstein, Ashley Chabot, Guolian Kang, and Shannon McKinney-Freeman.  Lifelong haematopoiesis is established by hundreds of precursors throughout mammalian ontogeny.  Nature cell biology, 19(10):1153–1163, 2017.
  * [11] Olga Gorlova, Bo Peng, David Yankelevitz, Claudia Henschke, and Marek Kimmel.  Estimating the growth rates of primary lung tumours from samples with missing measurements.  Statistics in medicine, 24(7):1117–1134, 2005.
  * [12] I. S. Gradshteyn and I. M. Ryzhik.  Table of integrals, series, and products.  Elsevier/Academic Press, Amsterdam, seventh edition, 2007.  Translated from the Russian, Translation edited and with a preface by Alan Jeffrey and Daniel Zwillinger, With one CD-ROM (Windows, Macintosh and UNIX).
  * [13] G. Grimmett and D. Stirzaker.  Probability and random processes.  Oxford university press, 2020.
  * [14] T. E. Harris.  Branching Processes.  The Annals of Mathematical Statistics, 19(4):474–494, 1948.
  * [15] Daniel Hormaechea-Agulla, Katie A Matatall, Duy T Le, Bailee Kain, Xiaochen Long, Pawel Kus, Roman Jaksik, Grant A Challen, Marek Kimmel, and Katherine Y King.  Chronic infection drives dnmt3a-loss-of-function clonal hematopoiesis via ifn$\gamma$ signaling.  Cell Stem Cell, 2021.
  * [16] Marek Kimmel.  Stochasticity and determinism in models of hematopoiesis.  In: A systems biology approach to blood; Corey, Kimmel and Leonard, eds., pages 119–152, 2014.
  * [17] Marek Kimmel and David Axelrod.  Branching Processes in Biology, volume 19.  Springer, New York, NY, 2002.
  * [18] Marek Kimmel and Seth Corey.  Stochastic hypothesis of transition from inborn neutropenia to aml: Interactions of cell population dynamics and population genetics.  Frontiers in Oncology, 3:89, 2013.
  * [19] C.A. Macken and A.S. Perelson.  Stem Cell Proliferation and Differentiation: A Multitype Branching Process Model.  Lecture Notes in Biomathematics. Springer Berlin Heidelberg, 1988.
  * [20] Michael C Mackey and Leon Glass.  Oscillation and chaos in physiological control systems.  Science, 197(4300):287–289, 1977.
  * [21] Patricia A Olofsen and Ivo P Touw.  Modeling severe congenital neutropenia in induced pluripotent stem cells.  In Recent Advances in iPSC Disease Modeling, Volume 1, pages 85–101. Elsevier, 2020.
  * [22] Laurent Pujo-Menjouet.  Blood cell dynamics: Half of a century of modelling.  Mathematical Modelling of Natural Phenomena, 11:92–115, 2016.
  * [23] Bryan E Shepherd, Peter Guttorp, Peter M Lansdorp, and Janis L Abkowitz.  Estimating human hematopoietic stem cell kinetics using granulocyte telomere lengths.  Experimental hematology, 32(11):1040–1050, 2004.
  * [24] Igor Sidorov, Masayuki Kimura, Anatoli Yashin, and Abraham Aviv.  Leukocyte telomere dynamics and human hematopoietic stem cell kinetics during somatic growth.  Experimental hematology, 37(4):514–524, 2009.
  * [25] T Stiehl, AD Ho, and A Marciniak-Czochra.  The impact of cd34+ cell dose on engraftment after scts: personalized estimates based on mathematical modeling.  Bone marrow transplantation, 49(1):30–37, 2014.
  * [26] Thomas Stiehl, Natalia Baran, Anthony D Ho, and Anna Marciniak-Czochra.  Clonal selection and therapy resistance in acute leukaemias: mathematical modelling explains different proliferation patterns at diagnosis and relapse.  Journal of The Royal Society Interface, 11(94):20140079, 2014.
  * [27] James E. Till, Ernest A. McCulloch, and Louis Siminovitch.  A stochastic model of stem cell proliferation, based on the growth of spleen colony-forming cells.  Proceedings of the National Academy of Sciences of the United States of America, 51:29–36, 1964.
  * [28] Ana Villalba, Magdalena Götz, and Víctor Borrell.  The regulation of cortical neurogenesis.  Current topics in developmental biology, 142:1–66, 2020.
  * [29] Caihong Wang, Whitney Trotter Ross, and Indira U Mysorekar.  Urothelial generation and regeneration in development, injury, and cancer.  Developmental Dynamics, 246(4):336–343, 2017.
  * [30] Zakary L. Whichard, Casim A. Sarkar, Marek Kimmel, and Seth J. Corey.  Hematopoiesis and its disorders: a systems biology approach.  Blood, 115(12):2339–2347, 2010.
  * [31] Tomasz Wojdyla, Hrishikesh Mehta, Taly Glaubach, Roberto Bertolusso, Marta Iwanaszko, Rosemary Braun, Seth J. Corey, and Marek Kimmel.  Mutation, drift and selection in single-driver hematologic malignancy: Example of secondary myelodysplastic syndrome following treatment of inherited neutropenia.  PLOS Computational Biology, 15(1):1–22, 2019.

Supplementary Materials for
“Stochastic Models of Stem Cells and Their Descendants under Different
Criticality Assumptions”

Nam H Nguyen and Marek Kimmel

Department of Statistics, Rice University, Houston, TX, USA

In this supplement, we will give detailed proofs that have been not covered in
the main article. We will frequently refer to equations and expressions in the
main text.

### Proof of Theorem 1b

Most of the calculations in the proof of Theorem 1a remain unchanged. We start
from the second-order ODE (4.2.10). Up to a multiplicative constant, the
general solution to this second-order linear ODE can be expresed in terms of
the Whittaker functions

$T(w)=M(\theta_{1},\theta_{2},\theta w)+CW(\theta_{1},\theta_{2},\theta w)$

where

$\theta_{1}=-\frac{\mu(1-\alpha)\alpha}{|1-2\alpha|};\
\theta_{2}=\frac{1}{2}\sqrt{1-4\mu\alpha+4\mu\alpha^{2}};\
\theta=|1-2\alpha|\mu$

Using formulae for the derivative of a Whittaker function (3), and short-hand
notations (1), we can show that

$T^{\prime}(w)=\Big{(}\frac{\theta}{2}-\frac{\theta_{1}}{w}\Big{)}(M_{0}(w)+CW_{0}(w))+\frac{1}{w}\Big{[}\Big{(}\frac{1}{2}+\theta_{1}+\theta_{2}\Big{)}M_{1}(w)-CW_{1}(w)\Big{]}$

We reverse the transformations using (4.2.11). Then, the explicit expression
of $F_{A}(w)$ is given by (1), where
$w=\frac{1}{1-f}=\frac{\lambda_{B}t(1-y)/4+1}{1-y}$ is a function of $t$.
Finally, using the initial condition $F_{A}\Big{(}w=\frac{1}{1-y}\Big{)}=x$,
we can deduce the following expression for the constant $C$.

$C=\frac{(\frac{1}{2}+\theta_{1}+\theta_{2}\big{)}M_{1}\big{(}\frac{1}{1-y}\big{)}-\kappa
M_{0}\big{(}\frac{1}{1-y}\big{)}}{\kappa
W_{0}\big{(}\frac{1}{1-y}\big{)}+W_{1}\big{(}\frac{1}{1-y}\big{)}}$

where

$\kappa=\frac{1}{1-y}\Big{(}-\mu(1-\alpha)^{2}x+\frac{\mu(1-2\alpha(1-\alpha)y)}{2}-\frac{\theta}{2}+\theta_{1}(1-y)\Big{)}$

### Proof of Theorem 1c

Combining (2.0.2) and (4.2.2), we obtain the following backward equation for
$F_{B}$

$\frac{dF_{B}}{dt}=\lambda_{B}(q^{2}F_{B}^{2}+2pqF_{B}-F_{B}+q^{2})$

Using basic integration, we get the following general solution

$F_{B}(x,y,t)=\frac{Pe^{\lambda_{B}(p-q)t}-p^{2}}{Pe^{\lambda_{B}(p-q)t}-q^{2}}$

Using the initial condition $F_{B}(x,y,0)=y$, we find that
$P=\frac{p^{2}-yq^{2}}{1-y}$. From (4.2.1), we need to solve the backward
equation for $F_{A}$

$\frac{dF_{A}}{dt}=-\lambda_{A}+\lambda_{A}((1-\alpha)F_{A}+\alpha F_{B})^{2}$

Write $z=Pe^{\lambda_{B}(p-q)t}$. Then, $F_{B}(z)=\frac{z-p^{2}}{z-q^{2}}$.
Note also that $\frac{dz}{dt}=\lambda_{B}(p-q)z$. Under this transformation,
we end up with the following equation for $F_{A}$ in terms of $z$

$\frac{dF_{A}}{dz}\lambda_{B}(p-q)z=-\lambda_{A}F_{A}+\lambda_{A}\bigg{(}(1-\alpha)F_{A}+\alpha\frac{z-p^{2}}{z-q^{2}}\bigg{)}^{2}$

Denoting $\mu=\frac{\lambda_{A}}{\lambda_{B}(p-q)}$,

$\frac{dF_{A}}{dz}=\frac{\mu}{z}(1-\alpha)^{2}F_{A}^{2}+\frac{\mu}{z}\bigg{(}2\alpha(1-\alpha)\frac{z-p^{2}}{z-q^{2}}-1\bigg{)}F_{A}+\frac{\mu}{z}\alpha^{2}\bigg{(}\frac{z-p^{2}}{z-q^{2}}\bigg{)}^{2}$

Let’s denote the coefficients of $F_{A}^{2}$, $F_{A}$ and the constant term by
$A$, $B$ and $C$ respectively

$\displaystyle A=\frac{\mu(1-\alpha)^{2}}{z};\
B=\frac{\mu}{z}\bigg{(}2\alpha(1-\alpha)\frac{z-p^{2}}{z-q^{2}}-1\bigg{)};\
C=\frac{\mu\alpha^{2}}{z}\bigg{(}\frac{z-p^{2}}{z-q^{2}}\bigg{)}^{2}$

As before, this non-linear ODE can be transformed to a Sturm-Liouville
differential equation

$u^{\prime\prime}+\gamma u^{\prime}+\beta u=0$

where $\gamma$ and $\beta$ are given by

$\gamma=-\bigg{(}\frac{A^{\prime}}{A}+B\bigg{)}=\frac{1}{z}-\frac{\mu}{z}\bigg{(}2\alpha(1-\alpha)\frac{z-p^{2}}{z-q^{2}}-1\bigg{)};\
\beta=AC=\bigg{(}\frac{\mu\alpha(1-\alpha)}{z}\bigg{)}^{2}\bigg{(}\frac{z-p^{2}}{z-q^{2}}\bigg{)}^{2}$

Then, we apply the transformation $u=\Phi T$, where $\Phi$ is a function of
$z$ that satisfies the condition $\Phi^{\prime}=-\gamma\Phi/2$. The equation
above reduces to

$T^{\prime\prime}+\bigg{[}\frac{1-\mu^{2}}{4z^{2}}+\frac{\mu\alpha(1-\alpha)}{z(z-q^{2})}\bigg{(}\frac{\mu(z-p^{2})}{z}+\frac{p-q}{z-q^{2}}\bigg{)}\bigg{]}T=0$

The solution to this ODE can be written in terms of the hypergeometric
functions

$T(z)=(z-q^{2})^{\frac{1+\theta_{2}}{2}}\bigg{(}z^{\frac{1+\theta_{1}}{2}}F_{0}^{+}(z)+Cz^{\frac{1-\theta_{1}}{2}}F_{0}^{-}(z)\bigg{)}$

where

$\theta_{1}=\frac{\mu\sqrt{q^{2}-4\alpha(1-\alpha)p^{2}}}{q};\
\theta_{2}=\frac{\sqrt{q^{2}+4\alpha(1-\alpha)\mu(q-p)}}{q};\
\theta_{3}=\mu|1-2\alpha|$

and $F_{0}^{+}$ and $F_{0}^{-}$ are short-hand notations defined in (1). Using
formulae for the first derivative of a hypergeometric function (3.0.8), we
have that

$\displaystyle T^{\prime}(z)=(z-q^{2})^{\frac{1+\theta_{2}}{2}}$
$\displaystyle\bigg{[}\frac{1+\theta_{1}}{2}z^{-\frac{1-\theta_{1}}{2}}F_{0}^{+}(z)+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}z^{\frac{1+\theta_{1}}{2}}F_{1}^{+}(z)$
$\displaystyle+C\bigg{(}\frac{1-\theta_{1}}{2}z^{-\frac{1+\theta_{1}}{2}}F_{0}^{-}(z)+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}z^{\frac{1-\theta_{1}}{2}}F_{1}^{-}(z)\bigg{)}\bigg{]}$
$\displaystyle+\frac{1+\theta_{2}}{2(z-q^{2})}T(z)$

By reversing the transformations using (4.2.11), we find that the explicit
expression of $F_{A}(z)$ is given by (1). The initial condition in terms of
$z$ is $F_{A}\Big{(}z=\frac{p^{2}-yq^{2}}{1-y}\Big{)}=x$. The constant $C$ is
given as follows.

$\displaystyle C=$
$\displaystyle\bigg{[}\bar{z}^{\frac{1+\theta_{1}}{2}}F_{0}^{+}(\bar{z})\kappa-\bigg{(}\frac{1+\theta_{1}}{2}\bar{z}^{-\frac{1-\theta_{1}}{2}}F_{0}^{+}(\bar{z})+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}\bar{z}^{\frac{1+\theta_{1}}{2}}F_{1}^{+}(\bar{z})\bigg{)}\bigg{]}$
$\displaystyle\bigg{/}\bigg{[}\bigg{(}\frac{1-\theta_{1}}{2}\bar{z}^{-\frac{1+\theta_{1}}{2}}F_{0}^{-}(\bar{z})+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}\bar{z}^{\frac{1-\theta_{1}}{2}}F_{1}^{-}(\bar{z})\bigg{)}-\bar{z}^{\frac{1-\theta_{1}}{2}}F_{0}^{-}(\bar{z})\kappa\bigg{]}$

where $\bar{z}=\frac{p^{2}-yq^{2}}{1-y}$, and

$\kappa=\frac{1-y}{p^{2}-yq^{2}}\bigg{(}\frac{-\mu(x+y)}{4}+\frac{1+\mu}{2}-\frac{(1+\theta_{2})(p^{2}-yq^{2})}{2(p-q)}\bigg{)}$

### Proof of Theorem 2b

The probability of extinction by time $t$, $E(t)$, is $F_{A}(0,0,t)$, where
$F_{A}$ is given by (1). From Theorem 1b, at $x=0$ and $y=0$,

$C_{0}\coloneqq
C|_{x=0,y=0}=\frac{(\frac{1}{2}+\theta_{1}+\theta_{2}\big{)}M_{1}(1)-\kappa_{0}M_{0}(1)}{\kappa_{0}W_{0}(1)+W_{1}(1)}$

where

$\kappa_{0}\coloneqq\kappa|_{x=0,y=0}=\frac{\mu}{2}-\frac{\theta}{2}+\theta_{1}$

and $\theta_{1}$, $\theta_{2}$ and $\theta$ are constants defined in terms of
the model parameters, $\mu$ and $\alpha$, only. Therefore, it follows that
$C_{0}$ is a finite constant.
To study the asymptotic properties of the system, we use the algebraic
relationship (3) between the Whittaker functions and the confluent
hypergeometric functions. Note that, in Theorem 1b, $w_{0}\coloneqq
w|_{x=0,y=0}=\frac{\lambda_{B}}{4}t+1$ at $x=0$ and $y=0$. Hence,
$w_{0}\to\infty$ as $t\to\infty$, and we can use the large-argument properties
(3) of ${\vphantom{F}}_{1}\kern-0.5ptF_{1}(.)$ and $U(.)$. Due to the presence
of the exponential term, it is clear that the
${\vphantom{F}}_{1}\kern-0.5ptF_{1}(.)$ dominates $U(.)$ in magnitude when the
argument is large. Thus, $M(.)$ dominates $W(.)$ in absolute value
asymptotically. This, coupled with the knowledge that $C_{0}$ is finite in
large time, leads to

$\frac{\big{(}\frac{1}{2}+\theta_{1}+\theta_{2}\big{)}M_{1}(w_{0})-C_{0}W_{1}(w_{0})}{M_{0}(w_{0})+C_{0}W_{0}(w_{0})}\sim\bigg{(}\frac{1}{2}+\theta_{1}+\theta_{2}\bigg{)}\frac{M_{1}(w_{0})}{M_{0}(w_{0})}$

Again, using (3), and the approximation of
${\vphantom{F}}_{1}\kern-0.5ptF_{1}(.)$ as given by (3),

$\frac{M_{1}(w_{0})}{M_{0}(w_{0})}\sim\frac{\theta_{2}-\theta_{1}-\frac{1}{2}}{\theta
w_{0}}$

Then, it follows that

$\displaystyle E(t)$
$\displaystyle=\frac{-1}{\mu(1-\alpha)^{2}}\bigg{(}\frac{\theta}{2}-\frac{\theta_{1}}{w_{0}}+\frac{1}{w_{0}}\frac{\big{(}\frac{1}{2}+\theta_{1}+\theta_{2}\big{)}M_{1}(w_{0})-C_{0}W_{1}(w_{0})}{M_{0}(w_{0})+C_{0}W_{0}(w_{0})}$
$\displaystyle\qquad\qquad\qquad-\frac{\mu(w_{0}+2\alpha(1-\alpha)(1-w_{0}))}{2w_{0}}\bigg{)}$
$\displaystyle\sim\frac{-1}{\mu(1-\alpha)^{2}}\bigg{(}\frac{\theta}{2}-\frac{\theta_{1}}{w_{0}}+\frac{1}{w_{0}}\frac{\big{(}\frac{1}{2}+\theta_{1}+\theta_{2}\big{)}\big{(}\theta_{2}-\theta_{1}-\frac{1}{2}\big{)}}{\theta
w_{0}}$
$\displaystyle\qquad\qquad\qquad-\frac{\mu(w_{0}+2\alpha(1-\alpha)(1-w_{0}))}{2w_{0}}\bigg{)}$
$\displaystyle\sim\frac{-1}{\mu(1-\alpha)^{2}}\bigg{(}\frac{\theta}{2}-\frac{\mu(1-2\alpha(1-\alpha))}{2}-\frac{\theta_{1}}{w_{0}}-\frac{\mu\alpha(1-\alpha)}{w_{0}}\bigg{)}$
$\displaystyle=\frac{-1}{\mu(1-\alpha)^{2}}\bigg{(}\frac{\mu|1-2\alpha|}{2}-\frac{\mu(1-2\alpha(1-\alpha))}{2}+\frac{\mu\alpha(1-\alpha)}{|1-2\alpha|w_{0}}-\frac{\mu\alpha(1-\alpha)}{w_{0}}\bigg{)}$
$\displaystyle=\begin{cases}\frac{\alpha^{2}}{(1-\alpha)^{2}}-\frac{2\alpha^{2}}{(1-\alpha)(1-2\alpha)w_{0}},&\mbox{if
}0<\alpha<\frac{1}{2}\\\ 1-\frac{2\alpha}{(2\alpha-1)w_{0}},&\mbox{if
}\frac{1}{2}<\alpha<1\end{cases}$
$\displaystyle\sim\begin{cases}\frac{\alpha^{2}}{(1-\alpha)^{2}}-\frac{8\alpha^{2}}{(1-\alpha)(1-2\alpha)\lambda_{B}t},&\mbox{if
}0<\alpha<\frac{1}{2}\\\ 1-\frac{8\alpha}{(2\alpha-1)\lambda_{B}t},&\mbox{if
}\frac{1}{2}<\alpha<1\end{cases}$

The theorem then follows.

### Proof of Theorem 2c

The probability of extinction by time $t$, $E(t)$, is $F_{A}(0,0,t)$, where
$F_{A}$ is given by (1). Recall from Theorem 1c that
$\theta_{1}=\frac{\mu\sqrt{q^{2}-4\alpha(1-\alpha)p^{2}}}{q}$. In the super-
critical case, we assume $p<q$ and, therefore, $\theta_{1}$ is real for all
$\alpha\in[0,1]$. Similarly, we can easily show that $\theta_{2}$ and
$\theta_{3}$ are real constants defined in terms of the model parameters for
all $\alpha$. From Theorem 1c, at $x=0$ and $y=0$,

$\bar{z}_{0}\coloneqq\bar{z}|_{x=0,y=0}=p^{2}$

where

$\kappa_{0}\coloneqq\kappa|_{x=0,y=0}=\frac{1}{p^{2}}\bigg{(}\frac{1+\mu}{2}-\frac{(1+\theta_{2}^{2})p^{2}}{2(p-q)}\bigg{)}$

are real and finite. It then follows that

$\displaystyle C_{0}$ $\displaystyle\coloneqq C|_{x=0,y=0}$
$\displaystyle=\bigg{[}\bar{z_{0}}^{\frac{1+\theta_{1}}{2}}F_{0}^{+}(\bar{z}_{0})\kappa_{0}-\bigg{(}\frac{1+\theta_{1}}{2}\bar{z}_{0}^{-\frac{1-\theta_{1}}{2}}F_{0}^{+}(\bar{z}_{0})+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}\bar{z}_{0}^{\frac{1+\theta_{1}}{2}}F_{1}^{+}(\bar{z}_{0})\bigg{)}\bigg{]}$
$\displaystyle\
\bigg{/}\bigg{[}\bigg{(}\frac{1-\theta_{1}}{2}\bar{z}_{0}^{-\frac{1+\theta_{1}}{2}}F_{0}^{-}(\bar{z}_{0})+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}\bar{z}_{0}^{\frac{1-\theta_{1}}{2}}F_{1}^{-}(\bar{z}_{0})\bigg{)}-\bar{z}_{0}^{\frac{1-\theta_{1}}{2}}F_{0}^{-}(\bar{z}_{0})\kappa_{0}\bigg{]}$

is finite. Recall that $z=\frac{p^{2}-yq^{2}}{1-y}e^{\lambda_{B}(p-q)t}$. At
$x=0$ and $y=0$, $z_{0}\coloneqq z|_{x=0,y=0}=p^{2}e^{\lambda_{B}(p-q)t}$.
Since $p<q$, we have that $z_{0}\to 0$ as $t\to\infty$. From the series
expansion of the hypergeometric function (3.0.7), it is easy to see that
hypergeometric functions are identically $1$ when evaluated at $0$. Setting
$w_{0}=\frac{1}{z_{0}}$, we have

$\displaystyle\bigg{[}\frac{1+\theta_{1}}{2}z_{0}^{-\frac{1-\theta_{1}}{2}}F_{0}^{+}(z_{0})+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}z_{0}^{\frac{1+\theta_{1}}{2}}F_{1}^{+}(z_{0})$
$\displaystyle+C_{0}\bigg{(}\frac{1-\theta_{1}}{2}z_{0}^{-\frac{1+\theta_{1}}{2}}F_{0}^{-}(z_{0})+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}z_{0}^{\frac{1-\theta_{1}}{2}}F_{1}^{-}(z_{0})\bigg{)}\bigg{]}$
$\displaystyle\bigg{/}\bigg{[}z_{0}^{\frac{1+\theta_{1}}{2}}F_{0}^{+}(z_{0})+C_{0}z_{0}^{\frac{1-\theta_{1}}{2}}F_{0}^{-}(z_{0})\bigg{]}$
$\displaystyle\sim\bigg{[}\frac{1+\theta_{1}}{2}w_{0}^{\frac{1-\theta_{1}}{2}}+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}w_{0}^{-\frac{1+\theta_{1}}{2}}$
$\displaystyle+C_{0}\bigg{(}\frac{1-\theta_{1}}{2}w_{0}^{\frac{1+\theta_{1}}{2}}+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}w_{0}^{-\frac{1-\theta_{1}}{2}}\bigg{)}\bigg{]}$
$\displaystyle\bigg{/}\bigg{[}w_{0}^{-\frac{1+\theta_{1}}{2}}+C_{0}w_{0}^{-\frac{1-\theta_{1}}{2}}\bigg{]}$

Clearly, $w_{0}\to\infty$ as $t\to\infty$, so we are interested in the
behavior of the expression above when $w_{0}$ is large. Since
$\mu=\frac{\lambda_{A}}{\lambda_{B}(p-q)}$ is negative when $p<q$, it follows
that $\theta_{1}$ is also negative. By comparing the exponents,
$w_{0}^{\frac{1-\theta_{1}}{2}}$ is the dominant term in the numerator, and
$w_{0}^{-\frac{1+\theta_{1}}{2}}$ is the dominant term in the denominator.
This, coupled with the fact that $C_{0}$ is finite, gives us the following
approximation

$\displaystyle\bigg{[}\frac{1+\theta_{1}}{2}w_{0}^{\frac{1-\theta_{1}}{2}}+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}w_{0}^{-\frac{1+\theta_{1}}{2}}$
$\displaystyle+C_{0}\bigg{(}\frac{1-\theta_{1}}{2}w_{0}^{\frac{1+\theta_{1}}{2}}+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}w_{0}^{-\frac{1-\theta_{1}}{2}}\bigg{)}\bigg{]}$
$\displaystyle\bigg{/}\bigg{[}w_{0}^{-\frac{1+\theta_{1}}{2}}+C_{0}w_{0}^{-\frac{1-\theta_{1}}{2}}\bigg{]}\sim\frac{1+\theta_{1}}{2}w_{0}$

Assuming $-1<\theta_{1}<0$, it follows that

$\displaystyle\bigg{[}\frac{1+\theta_{1}}{2}w_{0}^{\frac{1-\theta_{1}}{2}}+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}w_{0}^{-\frac{1+\theta_{1}}{2}}$
$\displaystyle+C_{0}\bigg{(}\frac{1-\theta_{1}}{2}w_{0}^{\frac{1+\theta_{1}}{2}}+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}w_{0}^{-\frac{1-\theta_{1}}{2}}\bigg{)}\bigg{]}$
$\displaystyle\bigg{/}\bigg{[}w_{0}^{-\frac{1+\theta_{1}}{2}}+C_{0}w_{0}^{-\frac{1-\theta_{1}}{2}}\bigg{]}\sim\frac{1+\theta_{1}}{2}w_{0}-C_{0}\theta_{1}w_{0}^{1+\theta_{1}}$

Then, for large $t$,

$\displaystyle E(t)$
$\displaystyle=\frac{-z_{0}}{\mu(1-\alpha)^{2}}\bigg{\\{}\bigg{[}\frac{1+\theta_{1}}{2}z_{0}^{-\frac{1-\theta_{1}}{2}}F_{0}^{+}(z_{0})+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}z_{0}^{\frac{1+\theta_{1}}{2}}F_{1}^{+}(z_{0})$
$\displaystyle+C_{0}\bigg{(}\frac{1-\theta_{1}}{2}z_{0}^{-\frac{1+\theta_{1}}{2}}F_{0}^{-}(z_{0})+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}z_{0}^{\frac{1-\theta_{1}}{2}}F_{1}^{-}(z_{0})\bigg{)}\bigg{]}$
$\displaystyle\bigg{/}\bigg{[}z_{0}^{\frac{1+\theta_{1}}{2}}F_{0}^{+}(z_{0})+C_{0}z_{0}^{\frac{1-\theta_{1}}{2}}F_{0}^{-}(z_{0})\bigg{]}+\frac{1}{2}\frac{1+\theta_{2}}{z_{0}-q^{2}}$
$\displaystyle-\frac{1}{2}\bigg{[}\frac{1}{z_{0}}-\frac{\mu}{z_{0}}\bigg{(}2\alpha(1-\alpha)\frac{z_{0}-p^{2}}{z_{0}-q^{2}}-1\bigg{)}\bigg{]}\bigg{\\}}$
$\displaystyle\sim\frac{-z_{0}}{\mu(1-\alpha)^{2}}\bigg{[}\frac{1+\theta_{1}}{2}\frac{1}{z_{0}}-C_{0}\theta_{1}\frac{1}{z_{0}^{1+\theta_{1}}}-\frac{1}{2}\frac{1+\theta_{2}}{q^{2}}-\frac{1}{2z_{0}}+\frac{\mu}{2z_{0}}\bigg{(}2\alpha(1-\alpha)\frac{p^{2}}{q^{2}}-1\bigg{)}\bigg{]}$
$\displaystyle\sim\frac{-z_{0}}{\mu(1-\alpha)^{2}}\bigg{[}\frac{\theta_{1}}{2z_{0}}+\frac{\mu}{2z_{0}}\bigg{(}2\alpha(1-\alpha)\frac{p^{2}}{q^{2}}-1\bigg{)}-C_{0}\theta_{1}\frac{1}{z_{0}^{1+\theta_{1}}}\bigg{]}$
$\displaystyle=\frac{1}{2(1-\alpha)^{2}}\bigg{(}1-\frac{2\alpha(1-\alpha)p^{2}}{q^{2}}-\frac{\sqrt{q^{2}-4\alpha(1-\alpha)p^{2}}}{q}\bigg{)}+\frac{C_{0}\theta_{1}p^{2}}{\mu(1-\alpha)^{2}}e^{-\lambda_{B}(p-q)\theta_{1}t}$

Again, it is straight-forward to show that the constant part is between $0$
and $1$ and, therefore, represents a valid probability of extinction.

The same approach works for the case $\theta_{1}<-1$. We can show that

$\displaystyle\bigg{[}\frac{1+\theta_{1}}{2}w_{0}^{\frac{1-\theta_{1}}{2}}+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}w_{0}^{-\frac{1+\theta_{1}}{2}}$
$\displaystyle+C_{0}\bigg{(}\frac{1-\theta_{1}}{2}w_{0}^{\frac{1+\theta_{1}}{2}}+\frac{(1-\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1-\theta_{1})}w_{0}^{-\frac{1-\theta_{1}}{2}}\bigg{)}\bigg{]}$
$\displaystyle\bigg{/}\bigg{[}w_{0}^{-\frac{1+\theta_{1}}{2}}+C_{0}w_{0}^{-\frac{1-\theta_{1}}{2}}\bigg{]}\sim\frac{1+\theta_{1}}{2}w_{0}+\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}$

Then, the probability of extinction in large time is approximated by

$\displaystyle E(t)$
$\displaystyle\sim\frac{1}{2(1-\alpha)^{2}}\bigg{(}1-\frac{2\alpha(1-\alpha)p^{2}}{q^{2}}-\frac{\sqrt{q^{2}-4\alpha(1-\alpha)p^{2}}}{q}\bigg{)}$
$\displaystyle-\frac{p^{2}}{\mu(1-\alpha)^{2}}\bigg{[}\frac{(1+\theta_{1}+\theta_{2})^{2}-\theta_{3}^{2}}{4q^{2}(1+\theta_{1})}-\frac{1+\theta_{2}}{2q^{2}}\bigg{]}e^{-\lambda_{B}(p-q)\theta_{1}t}$

The theorem then follows.