# Immortal Branching Processes

P. L. Krapivsky<EMAIL_ADDRESS>Department of Physics, 590 Commonwealth Avenue,
Boston University, Boston, MA 02215, USA  S. Redner<EMAIL_ADDRESS>Santa
Fe Institute, 1399 Hyde Park Road, Santa Fe, New Mexico 87501, USA

###### Abstract

We investigate the dynamics of an _immortal_ branching process. In the
classic, critical branching process, particles give birth to a single
offspring or die at the same rates. Even though the average population is
constant in time, the ultimate fate of the population is extinction. We
augment this branching process with immortality by positing that either: (a) a
single particle cannot die, or (b) there exists an immortal stem cell that
gives birth to ordinary cells that can subsequently undergo critical
branching. We also treat immortal two-type branching. We discuss the new
dynamical aspects of these immortal branching processes.

††journal: Physica A

##  1 Introduction

We dedicate this contribution to the memory of Dietrich Stauffer, a singular
personality in the physics community. He was a man with strong views, who
spoke with few compunctions, and who had a great sense of humor. Looking back
at the emails that he exchanged with one of us (SR) brought back pleasant
memories of his caustic wit. One such example is Dietrich’s email to another
person (with a cc to me) in which he writes: “I am very grateful that you try
to educate the Canadian citizen Sid Redner [yes, SR is Canadian] in a
civilised language. But since he is a traitor working for the Southern Enemy
[meaning the USA], I doubt that you will have a long-lasting success. He does
not even use the spelling of Her Britannic Majesty, Canada’s head of state.”
He enjoyed making fun of SR’s allegedly traitorous behavior (or behaviour) and
SR enjoyed getting poked by him. In addition to the many verbal and written
barbs that we exchanged, it was also a great pleasure to discuss science and
indeed any field of academic study with him. He will be sorely missed.

Dietrich is best known for his contributions to the field of percolation [1].
Percolation has applications to diverse fields; one such example is epidemic
spread. From the dynamical perspective, the temporal course of an epidemic can
be viewed as a birth/death or a branching process. Namely, an infected
individual can infect someone who is not yet sick and an infected individual
can be “removed” by either recovering or expiring. If we denote an infected
individual as P, the epidemic dynamics can be symbolically represented as

$P\to P+P$ | rate $~{}r$
---|---
$P\to\emptyset$ | rate $~{}1$
(1)

where the death rate is set to one by choosing appropriate time units. The
_critical_ branching process in which the birth and death rates are the same,
$r=1$, is especially popular, and has numerous applications [2, 3, 4].
Although the adjective “critical” suggests that this branching process is
idiosyncratic because of the tuning of the birth and death rates, this case is
arguably more important than its subcritical and supercritical brethren.
Indeed, subcritical branching ($r<1$) quickly goes extinct while supercritical
branching ($r>1$) results in unlimited growth.

Branching processes are related to the stochastic version of the susceptible-
infected-recovered (SIR) infection model, in which Ps correspond to infected
individuals and $\emptyset$s to recovered. In the SIR model [5, 6, 7, 8], the
population consists of susceptible, infected, and recovered individuals. The
infection spreads by contact between infected and susceptible individuals,
while infected individuals spontaneously recover or die. In the context of
infectious diseases, critical branching processes arise as a balance between
human efforts that strive to reduce the infection rate [6] and natural
evolution that increases the infectivity of diseases [7]. Critical branching
can also describe homeostasis (the maintenance of healthy tissues) [9, 10].
For instance, the maintenance of skin cells arises when progenitor cells (P),
which exist far below the epidermis, undergo critical branching and (instead
of death) differentiate into post-mitotic cells that are eventually shed.

In this contribution, we study _immortal_ branching processes (IBPs). Our
first example is identical to the critical branching process with a simple
twist—when a single cell remains, it cannot die. One of the motivations for
this type of modeling is the long-term dynamics of HIV infection in which a
nearly undetectable level of viral load can linger in patients for a decade or
longer and then subsequently increase [11]. This branching without the
possibility of extinction greatly oversimplifies the HIV dynamics, but perhaps
it points the way to more realistic models. We analyze this model in Sec. 3
and demonstrate unexpectedly subtle dynamics that are replete with logarithmic
corrections in the long time limit.

Another IBP is a critical branching process that is augmented by the presence
of a single immortal stem cell whose properties are not affected by progenitor
cells. The reaction steps of this IBP are

$S\to S+P$ | rate $~{}\beta$
---|---
$P\to P+P$ | rate $~{}1$
$P\to\emptyset$ | rate $~{}1$
(2)

The first process accounts for a stem cell (S) that creates a single
progenitor cell P, while the latter two reactions represent critical
branching. Mathematically, the reactions (2) correspond to critical branching
with input. The strength of the source, which is quantified by $\beta$, can
vary with time in principle. For simplicity, we focus on the case where
$\beta=\text{const}$. There is no extinction because the stem cell is immortal
and the average number of progenitor cells grows as $\beta t$.

##  2 Critical Branching

It is helpful to first recall basic facts about critical branching. The
probability $P_{m}(t)$ to have $m$ cells at time $t$ satisfies the rate
equation

$\frac{dP_{m}}{dt}=(m+1)P_{m+1}-2mP_{m}+(m-1)P_{m-1}\,.$ (3)

The first term on the right-hand side corresponds to any one of $m+1$ cells
dying with rate 1; the other terms have similar explanations. We always assume
that the system initially contains a single cell, $P_{m}(t=0)=\delta_{m,1}$.

There are various ways to solve Eqs. (3). Perhaps the simplest makes use of
the exponential ansatz (see, e.g., [12])

$P_{m}(t)=A(t)[a(t)]^{m-1},$ (4)

which often works for infinite sets of equations with linear dependences on
$m$. With ansatz (4), the infinite set of equations (3) reduces to
$\frac{dA}{dt}=2A(a-1)$ and $\frac{da}{dt}=(a-1)^{2}$, from which
$A=(1+t)^{-2}$ and $a=t/(1+t)$. Thus

$P_{m}(t)=\frac{t^{m-1}}{(1+t)^{m+1}}$ (5)

is the probability that there are $m$ cells at time $t$. To appreciate the
meaning of this solution, it is useful to look at the integer moments

$\langle m^{k}\rangle=\sum_{m\geq 1}m^{k}P_{m}(t)\,.$

The first few moments are $\langle m^{0}\rangle=(1+t)^{-1},~{}\langle
m^{1}\rangle=1,~{}\langle m^{2}\rangle=1+2t$. The zeroth moment is just the
probability that the population does not go extinct—the survival probability.
Even though $\langle m^{0}\rangle\to 0$ as $t\to\infty$ so that eventual
extinction necessarily occurs, the mean extinction time is infinite. The
conserved nature of the first moment reflects the equality of the birth and
death rates so that the average population size is constant. In contrast, for
those realizations that do not go extinct, the population grows linearly with
time, as manifested by the behavior of the second moment.

##  3 Branching Without Extinction

We now amend critical branching so that extinction is impossible. A simple way
to impose this constraint is to retain the rules of critical branching but
forbid death when only a single cell remains. That is, Eqs. (3) remain valid
for $m\geq 2$, while the probability $P_{1}(t)$ now evolves according to

$\frac{dP_{1}}{dt}=-P_{1}+2P_{2}\,.$ (6)

To solve Eqs. (3) and (6), first notice that the ansatz (4) no longer solves
them, which suggests that the solution does not have an exponential form.
Instead, we proceed conventionally and introduce the Laplace transform

$\widetilde{P}_{m}(s)=\int_{0}^{\infty}dt\,e^{-st}P_{m}(t)\,,$

to recast (3) and (6) into the linear algebraic equations

$\displaystyle\begin{split}(s+2m)\widetilde{P}_{m}&=(m-1)\widetilde{P}_{m-1}+(m+1)\widetilde{P}_{m+1}\qquad
m\geq 2\,,\\\ (s+1)\widetilde{P}_{1}&=2\widetilde{P}_{2}+1\,.\end{split}$ (7)

We now apply the generating function method. By standard steps, we find that
the generating function

$\mathcal{P}(s,z)=\sum_{m\geq 1}\widetilde{P}_{m}(s)\,z^{m}$

satisfies the ordinary differential equation (ODE)

$s\mathcal{P}=(z-1)^{2}\,\frac{d\mathcal{P}}{dz}+(z-1)\widetilde{P}_{1}(s)-z\,.$
(8)

The details of the solution to this equation and the inversion of the Laplace
transform are given in A. The final result for the population distribution in
the Laplace domain is

$\widetilde{P}_{m}(s)=\frac{1}{se^{s}\Gamma(0,s)}\int_{0}^{\infty}d\eta\,e^{-s\eta}\,\frac{\eta^{m-1}}{(1+\eta)^{m+1}}\,.$
(9)

While it does not seem possible to invert this Laplace transform and obtain
$P_{m}(t)$ in terms of elementary or standard special functions, we can
determine the long-time behavior from the $s\to 0$ asymptotic of the Laplace
transform. Using the asymptotic formula (see, e.g., [13])

$\Gamma(0,s)=-\ln s-\gamma+\sum_{k\geq 1}\frac{(-s)^{k}}{k\times k!}\,,$

where $\gamma=0.577\,215\ldots$ is Euler’s constant and
$\Gamma(0,s)=\int_{s}^{\infty}du\,e^{-u}/u$ is the incomplete Gamma function,
the leading $s\to 0$ behavior of $\widetilde{P}_{1}(s)$ is

$\widetilde{P}_{1}(s)\simeq\frac{1}{s(-\ln s-\gamma)}\,,$ (10a) from which we
deduce the unexpectedly slow long-time decay $P_{1}(t)\simeq\frac{1}{\ln
t}\,.$ (10b)

This behavior strongly contrasts with $P_{1}(t)\simeq t^{-2}$ in critical
branching. Because there is no mechanism for an isolated particle to die, this
leads to a much slower temporal decay of $P_{1}(t)$.

We now use this solution for $P_{1}(t)$ in conjunction with the rate equations
(3) and (6) to recursively find the asymptotic behavior of $P_{m}(t)$.
Substituting (10b) into (6) leads to

$P_{2}(t)\simeq\frac{1}{2}\,\frac{1}{\ln t}$

to leading order. Substituting these asymptotic forms for $P_{1}(t)$ and
$P_{2}(t)$ into the rate equation (3) with $m=2$, we find

$P_{3}(t)\simeq\frac{1}{3}\,\frac{1}{\ln t}$

to leading order. Continuing this procedure iteratively gives

$P_{m}(t)\simeq\frac{1}{m}\,\frac{1}{\ln t}\,,$ (11)

which should be valid for $t\gg 1$ and $m\ll t$.

To determine the distribution $P_{m}(t)$ when $m\sim t$, it is helpful to
first determine the asymptotic growth of the moments $\langle m^{k}\rangle$.
Using Eqs. (3) and (6), we deduce the rate equation $\frac{d\langle
m\rangle}{dt}=P_{1}$. Combining this with the expression for $P_{1}(t)$ in
(10b), the asymptotic growth of the average number of cells is

$\langle m\rangle\simeq\frac{t}{\ln t}\,.$ (12)

Similarly $\frac{d\langle m^{2}\rangle}{dt}=2\langle m\rangle+P_{1}$ which in
conjunction with (10b) and (12) gives

$\langle m^{2}\rangle\simeq\frac{t^{2}}{\ln t}\,.$ (13)

For general $k\geq 2$, the evolution equation for $\langle m^{k}\rangle$ is

$\frac{d\langle
m^{k}\rangle}{dt}=2\sum_{a=1}^{\left\lfloor\frac{k}{2}\right\rfloor}\binom{k}{2a}\langle
m^{k-2a+1}\rangle+P_{1}\,,$ (14a) whose leading behavior is $\frac{d\langle
m^{k}\rangle}{dt}\simeq 2\binom{k}{2}\langle m^{k-1}\rangle\,.$ (14b)

The results for the first two moments suggest that $\langle m^{k}\rangle\simeq
C_{k}\,t^{k}/\ln t$. Substituting this ansatz into (14b), we find that this
ansatz is consistent and moreover it serves to fix the amplitude $C_{k}$. The
final result is

$\langle m^{k}\rangle\simeq(k-1)!\;\frac{t^{k}}{\ln t}$ (15)

for $k\geq 1$. (The zeroth moment is $\langle m^{0}\rangle=\sum_{m\geq
1}P_{m}(t)=1$ because extinction is impossible.)

The small-$m$ tail of $P_{m}(t)$ and the behavior of the moments suggest that
the full distribution $P_{m}(t)$ has the following scaling behavior

$P_{m}(t)\simeq\frac{1}{m}\;\frac{1}{\ln t}\;F(\mu),\qquad\mu=\frac{m}{t}\,,$
(16)

with $F(0)=1$. Using this distribution, we now express the moments $\langle
m^{k}\rangle$ through the integrals of the scaled distribution:

$\langle m^{k}\rangle\simeq\frac{t^{k}}{\ln
t}\,\int_{0}^{\infty}d\mu\,\mu^{k-1}F(\mu)\,.$

Thus we have the condition $\int_{0}^{\infty}d\mu\,\mu^{k-1}F(\mu)=(k-1)!$,
which allows us to fix the scaling function: $F(\mu)=e^{-\mu}$. Finally, we
obtain

$P_{m}(t)\simeq\frac{1}{m}\,\frac{1}{\ln t}\,e^{-m/t}\,.$ (17)

It is instructive to compare this scaling form with that of the critical
branching process. From (5), this latter scaling form is

$\displaystyle P_{m}(t)\simeq\frac{1}{t^{2}}\,e^{-m/t}\,.$ (18)

We again see that the population distribution has a much slower temporal decay
than in critical branching; this occurs because of the effective “source” at
$m=1$.

##  4 Branching With Input

A biologically motivated version of immortal branching is based on stem cells
that are always active. We focus on a system that begins with a single stem
cell S. This stem cell can give birth to ordinary (mortal) cells,
$\text{S}\to\text{S}+\text{P}$, and these mortal cells subsequently undergo
critical branching. We denote the rate at which stem cells give birth as
$\beta$; for ordinary cells, the birth and death rates are set to 1.

Denote again by $P_{m}(t)$ the probability to have $m$ cells, i.e., the stem
cell and $m-1$ mortal cells. The probability that only a stem cell exists
obeys

$\frac{dP_{1}}{dt}=-\beta P_{1}+P_{2}\,,$ (19a) while for $m\geq 2$ the rate
equations are $\displaystyle\frac{dP_{m}}{dt}$
$\displaystyle=mP_{m+1}-[2(m-1)+\beta]P_{m}+(m-2+\beta)P_{m-1}\,.$ (19b)

One can verify by direct substitution that

$P_{m}(t)=\frac{\Gamma(m-1+\beta)}{\Gamma(m)\,\Gamma(\beta)}\;\frac{t^{m-1}}{(1+t)^{m-1+\beta}}$
(20)

is the solution of Eqs. (19). In the scaling region $m,t\to\infty$ with
$\mu={m}/{t}$ finite, this solution acquires the scaling form

$P_{m}(t)\simeq
t^{-1}\Phi(\mu),\qquad\Phi(\mu)=\frac{\mu^{\beta-1}}{\Gamma(\beta)}\,e^{-\mu}\,.$
(21)

With this scaling form, the leading behavior of the moments are given by

$\langle m^{k}\rangle\simeq t^{k}\,\frac{\Gamma(\beta+k)}{\Gamma(\beta)}\,.$
(22)

The IBP described by Eqs. (19) is often known as a branching process with
immigration. The solution (20) was derived by Kendall [14]; see [15] for
recent work. We now discuss a multistage IBP that may have relevance to
cancer, where tumor cells often undergo multiple stages [16, 17]. The simplest
step in this direction to allow for the possibility of two-type branching
[18]. These models are analytically tractable [19]; models with three-type
branching have not been solved so far.

Two-type branching processes also give a good description of experimental data
[9] on the cell dynamics of skin tissue. The experimental data are well fit by
a model that involves progenitor cells (P) that can divide (proliferate) and
differentiate, and post-mitotic cells (M) that eventually disappear (leave the
basal layer). The elemental steps of this process are

$S\to S+P$ |  rate $\beta$
---|---
$P\to P+P$ |  rate $r$
$P\to P+M$ |  rate $1-2r$
$P\to M+M$ |  rate $r$
$M\to\emptyset$ |  rate $\gamma$
(23)

With no source ($\beta=0$), this two-type branching process was solved in
[20]. We now generalize the approach of Ref. [20] to $\beta>0$. We account for
the populations of progenitor and post-mitotic cells by $P_{m,n}(t)$, the
probability to have $m$ clones of $P$ and $n$ clones of $M$. This probability
distribution evolves according to the master equation

$\displaystyle\frac{dP_{m,n}}{dt}$
$\displaystyle=[r(m-1)+\beta]P_{m-1,n}+r(m+1)P_{m+1,n-2}$
$\displaystyle+(1-2r)mP_{m,n-1}+\gamma(n+1)P_{m,n+1}-(m+\gamma
n+\beta)P_{m,n}\,.$ (24)

We again introduce the generating function

$\mathcal{P}(x,y,t)=\sum_{m,n=0}^{\infty}x^{m}y^{n}P_{m,n}(t)\,,$

which satisfies the partial differential equation (PDE)

$\partial_{t}\mathcal{P}+U\partial_{x}\mathcal{P}+V\partial_{y}\mathcal{P}=W\mathcal{P}$
(25)

where $U=x(1-y)-r(x-y)^{2},~{}V=\gamma(y-1)$ and $W=\beta(x-1)$. For an
initially empty system $\mathcal{P}(x,y,t=0)=1$. Equation (25) is a first-
order hyperbolic PDE and like all such equations, can be solved using the
method of characteristics [21]. The characteristics are determined by
equations

$\displaystyle\begin{split}\frac{dx}{dt}&=U=x(1-y)-r(x-y)^{2}\\\
\frac{dy}{dt}&=V=\gamma(y-1)\end{split}$ (26)

By construction of the characteristic equations, the time derivative along a
characteristic is given by

$\frac{d\mathcal{P}}{dt}=\partial_{t}\mathcal{P}+U\partial_{x}\mathcal{P}+V\partial_{y}\mathcal{P}$

Combining this with Eq. (25) we get $\frac{d\mathcal{P}}{dt}=W\mathcal{P}$,
which we integrate to yield

$\mathcal{P}(X,Y,T)=\exp\\!\Big{\\{}\beta\int_{0}^{T}d\tau\,[x(t)-1]\Big{\\}}\,.$
(27)

It is convenient to think of $t$ as a ‘running’ time, so we use the notation
$T$ for the final time. The functions $x(t)$ and $y(t)$ vary along the
characteristic and rather than parameterizing by initial values $x(0)$ and
$y(0)$ it is preferable to parameterize by the final values $X=x(T)$ and
$Y=y(T)$.

To compute the generating function, Eq. (27), we need to determine $x(t)$
along the characteristic. We integrate the second of Eqs. (26) to give

$y-1=(Y-1)e^{\gamma(t-T)}\,.$ (28)

Substituting this into the first of (26) and changing the time variable $t$ to

$u=\gamma(1-Y)\,e^{\gamma(t-T)}\,,$ (29)

we arrive at

$\frac{dx}{du}=\frac{x}{\gamma}-\frac{r}{\gamma u}(x-1+u)^{2}\,.$ (30)

This Riccati equation admits an exact solution in terms of confluent
hypergeometric functions [20]. Thus one can obtain a solution, but it is
complicated and difficult to use to obtain explicit results in the general
case. A solution in terms of elementary functions is possible in the special
case $r=\frac{1}{4}$ and $\gamma=1$ (B). Two prominent results from this
analysis are: (i) $P_{m}(t)$, the probability to have $m$ progenitor cells,
and any number of post-mitotic cells and (ii) $P_{m,0}(t)$, the probability to
have $m$ progenitor cells, and _no_ post-mitotic cells. The former is given by

$P_{m}(t)=\frac{\Gamma(m+4\beta)}{\Gamma(4\beta)\,\Gamma(m+1)}\,\frac{1}{(1+t/4)^{4\beta}}\,\left(\frac{t/4}{1+t/4}\right)^{m}\,.$
(31a) This result coincides with the exact solution (20) with the
understanding that $m$ counts the number of progenitor cells, and the birth
and death rates are $r=1/4$. The latter probability distribution is
$P_{m,0}(t)=\frac{\Gamma(m+4\beta)}{\Gamma(4\beta)\,\Gamma(m+1)}\,\frac{e^{\beta(1-e^{-t})}}{(1+t/2)^{4\beta}}\,\left(\frac{t/4}{1+t/2}\right)^{m}\,.$
(31b)

##  5 Discussion

We analyzed two immortal branching processes (IBPs). The first is a simple
extension of critical branching in which extinction cannot occur. The emergent
behaviors are remarkably subtle, that are replete with logarithmic
corrections. The second is a branching process with a steady input of cells,
or equivalently, immigration. Such a steady input arises in a variety of non-
equilibrium many-body processes, such as aggregation with a steady input of
monomers [12], fragmentation with a steady input of large clusters [22], and
turbulence with a steady energy input at large scales [23, 24]. In these
examples, many new features were uncovered by incorporating steady input. Rich
behaviors also occur in spatially extended systems with a spatially localized
input [25, 26] and the role of spatial degrees of freedom in IBPs poses
interesting challenges.

## Acknowledgments

We benefited from correspondence with Tibor Antal. PLK thanks the hospitality
of the Santa Fe Institute, where this work began. The research of SR is
partially supported by NSF grant DMR-1910736.

##  Appendix A Derivation of Eq. (9)

To solve the governing differential equation (8) for the generating function,
we introduce the auxiliary variable $\zeta=1/(1-z)$ to recast this equation as

$s\mathcal{P}=\frac{d\mathcal{P}}{d\zeta}+1-\frac{\widetilde{P}_{1}(s)+1}{\zeta}\,.$
(32)

This linear inhomogeneous ordinary differential equation may be readily solved
to yield

$\mathcal{P}=\frac{1}{s}-\big{(}\widetilde{P}_{1}(s)+1\big{)}\int_{0}^{\infty}\frac{d\eta}{\eta+\zeta}\,e^{-s\eta}+\Phi(s)\,e^{s\zeta}\,.$
(33)

To fix the integration constant $\Phi(s)$ we take the Laplace transform of the
normalization condition $\sum_{m\geq 1}P_{m}(t)=1$ to give $\sum_{m\geq
1}\widetilde{P}_{m}(s)=1/s$. This then implies that $\lim_{z\to
1^{-}}\mathcal{P}(s,z)=1/s$. Since the limit $z\to 1^{-}$ corresponds to
$\zeta\to+\infty$, the second term on the right-hand side of (33) vanishes
while the third term, $\Phi(s)\,e^{s\zeta}$, diverges if $\Phi(s)\neq 0$. We
thus conclude that $\Phi(s)=0$.

Returning to the original variable $z$, we rewrite (33) as

$\mathcal{P}=\frac{1}{s}-\big{[}\widetilde{P}_{1}(s)+1\big{]}\int_{0}^{\infty}d\eta\;e^{-s\eta}\,\frac{1-z}{1+\eta-\eta
z}\,.$ (34)

We now expand in a Taylor series in $z$ to give the population distribution

$\widetilde{P}_{m}(s)=\big{[}\widetilde{P}_{1}(s)+1\big{]}\int_{0}^{\infty}d\eta\;e^{-s\eta}\,\frac{\eta^{m-1}}{(1+\eta)^{m+1}}\,.$
(35)

Specializing (35) to $m=1$ we obtain a closed equation for
$\widetilde{P}_{1}(s)$. After some straightforward steps, this equation
simplifies to

$\widetilde{P}_{1}(s)=\frac{1}{se^{s}\Gamma(0,s)}-1\,,$ (36)

where $\Gamma(0,s)=\int_{s}^{\infty}du\,e^{-u}/u$ is the incomplete Gamma
function. Substituting the above form for $\widetilde{P}_{1}(s)$ in (35), we
obtain Eq. (9).

##  Appendix B Solution for $P_{m}(t)$ and $P_{m,0}(t)$ when $r=\frac{1}{4}$
and $\gamma=1$

For $r=\frac{1}{4}$ and $\gamma=1$, the Riccati equation (30) admits a
solution [20] in terms of elementary functions

$x=1+u+\frac{1}{\frac{1}{4}\ln u+C}\,.$ (37)

When $t=T$, we have $x=X$ and $u=1-Y$ for $\gamma=1$. These requirements fix
the constant $C$ to be

$C=-\tfrac{1}{4}\,\ln(1-Y)-(2-X-Y)^{-1}\,.$ (38)

Using this result we rewrite (37) as

$x(t)=1+(1-Y)\,e^{t-T}-\left(\frac{T-t}{4}+\frac{1}{2-X-Y}\right)^{-1}\,.$
(39)

Substituting (39) into (27), computing the integral, and returning to the
original variables $x,y,t$, we obtain

$\mathcal{P}(x,y,t)=\frac{\exp\\!\big{[}\beta(1-e^{-t})(1-y)\big{]}}{\left[1+\frac{t}{4}(2-x-y)\right]^{4\beta}}\,.$
(40)

We now specialize (40) to the case $y=1$ and expand this generating function
in a Taylor series in $x$. This gives Eq. (31a). In a complementary way, we
specialize (40) to the case $x=1$ and use the definition

$\mathcal{P}(1,y,t)=\sum_{n=0}^{\infty}y^{n}\Pi_{n}(t)\,,$ (41)

which expresses this restricted generating function in terms of $\Pi_{n}(t)$,
the probability to have $n$ post-mitotic cells. Expanding $\mathcal{P}(1,y,t)$
in a Taylor series in $y$ leads to

$\Pi_{n}=e^{\beta(1-e^{-t})}\sum_{k=0}^{n}\frac{[-\beta(1-e^{-t})]^{k}}{k!}\,P_{n-k}(t)\,,$
(42)

with $P_{m}$ given by (31a).

Using (40) we obtain a simple expression for the probability to have no cells
of any kind:

$P_{0,0}(t)=\mathcal{P}(0,0,t)=\frac{e^{\beta(1-e^{-t})}}{(1+t/2)^{4\beta}}\,.$
(43)

The temporal behavior of (43) is the same as that for the probability to have
no progenitor cells, $P_{0}=(1+t/4)^{-4\beta}$ (Eq. (31a)), although with
different amplitudes: $P_{0}>P_{0,0}(t)$. Finally, we specialize (40) to $y=0$
to give

$\mathcal{P}(x,0,t)=\sum_{m=0}^{\infty}x^{m}P_{m,0}(t)=\frac{\exp\\!\big{[}\beta(1-e^{-t})\big{]}}{\left[1+\frac{t}{4}(2-x)\right]^{4\beta}}\,.$
(44)

Expanding this expression in a Taylor series gives the probabilities of having
only progenitor cells quoted in Eq. (31b).

## References

  * [1] D. Stauffer and A. Aharony, _Introduction to Percolation Theory_ , 2nd ed. (CRC press, 1994).
  * [2] W. Feller, An Introduction to Probability Theory and Its Applications, Vol. I, 3rd ed. (John Wiley, New York, 1968).
  * [3] T. E. Harris, The Theory of Branching Processes (Dover, New York, 1989).
  * [4] P. Haccou, P. Jagers, and V. A. Vatutin, Branching processes: variation, growth, and extinction of populations (Cambridge University Press, New York, 2005).
  * [5] N. T. J. Bailey, The Mathematical Theory of Infectious Diseases (Oxford University Press, Oxford, 1987).
  * [6] R. Anderson and R. May, Infectious Diseases: Dynamics and Control (Oxford University Press, Oxford, 1991).
  * [7] R. Antia, R. R. Regoes, J. C. Koella, and C. T. Bergstrom, Nature 426, 658 (2003).
  * [8] T. Lux, Physica A 567, 125710 (2021).
  * [9] E. Clayton, D. P. Doupé, A. M. Klein, D. J. Winton, B. D. Simons, and P. H. Jones, Nature 446, 185 (2007).
  * [10] P. Jones and B. D. Simons, Nat. Rev. Mol. Cell Biol. 9, 83 (2008).
  * [11] X. Wei et al., Nature 373, 117 (1995); M. A. Stafford et al., J. Theor. Biol. 203 295 (2000).
  * [12] P. L. Krapivsky, S. Redner and E. Ben-Naim, _A Kinetic View of Statistical Physics_ (Cambridge University Press, Cambridge, UK, 2010).
  * [13] C. M. Bender and S. A. Orszag, _Advanced Mathematical Methods for Scientists and Engineers_ (McGraw-Hill, New York, 1978).
  * [14] D. G. Kendall, Biometrika 35, 6 (1948).
  * [15] M. D. Nicholson and T. Antal, Bull. Math. Biol. 78, 2243 (2016).
  * [16] P. Armitage and R. Doll, Brit. J. Cancer 8, 1 (1954).
  * [17] N. Beerenwinkel et al., PLOS Comput. Biol. 3, e225 (2007).
  * [18] D. G. Kendall, Biometrica 47, 13 (1960).
  * [19] T. Antal and P. L. Krapivsky, J. Stat. Mech. P08018 (2011).
  * [20] T. Antal and P. L. Krapivsky, J. Stat. Mech. P07028 (2010).
  * [21] J. D. Logan, An Introduction to Nonlinear Partial Differential Equations (John Wiley& Sons, Inc., New Jersey, 2008).
  * [22] E. Ben-Naim and P.L. Krapivsky, Phys. Lett. A 275, 48 (2000).
  * [23] G. I. Barenblatt, _Scaling, Self-similarity, and Intermediate Asymptotics_ , (Cambridge University Press, Cambridge, UK, 1996).
  * [24] N. Goldenfeld and H. Y. Shih, J. Stat. Phys. 167, 575 (2017).
  * [25] Z. Cheng, S. Redner, and F. Leyvraz, Phys. Rev. Lett. 62, 2321 (1989).
  * [26] P. L. Krapivsky, J. Phys. A 48, 205003 (2015); J. Phys. A 48, 245002 (2015).