2. The Crump–Mode–Jagers branching process of clusters

Although we introduced the epidemic model from the perspective of individuals, it will be convenient for its analysis to look at clusters and their evolutions as time passes. Specifically, imagine that we create an (unoriented) edge between the infector and the infected at the time when a contamination occurs; each edge is further labeled traceable or untraceable, depending on the type of contamination. If we ignore the labels of edges, this endows the infected population at any given time with a genealogical tree structure which is rooted at the ancestor. Plainly this tree structure grows as new individuals are contaminated and new traceable or untraceable edges are added. The reader may find Figure 1 useful to visualize the notions that will be introduced.

Figure 1. Graphical representation of the epidemic at a given time. The ancestor is the vertex at the bottom of the figure. Vertices in red represent contagious individuals, vertices in white individuals who have been isolated. Full edges indicate traceable contaminations, and dotted ones untraceable contaminations. Clusters consist of subsets of vertices connected by full edges. In turn, clusters are connected by dotted edges. There are four active clusters: two at the first generation with sizes 4 and 3, one at the third generation with size 1, and one at the fourth generation with size 1. There are four isolated clusters: the ancestor cluster with size 5, one cluster with size 2 at the first generation, and two clusters with sizes 4 and 3 at the second generation.

Any pair of infected individuals is connected by a unique segment in the tree, which we call the contamination path. Two individuals belong to the same cluster if and only if their contamination path contains only traceable edges, and more generally, the number of untraceable edges along the contamination path between two individuals only depends on the two clusters to which these individuals belong. We then define the generation of a given cluster as the number of untraceable edges along the contamination path between any individual in that cluster and the ancestor.

We next observe that backtracking contaminations endows clusters with a natural genealogy, which in turn enables us to view the epidemic model as a so-called Crump–Mode–Jagers branching process; see [Reference Jagers9, Chapter 6] as well as [Reference Jagers10], [Reference Jagers and Nerman11], and [Reference Nerman16] for classical background, and also [Reference Holmgren and Janson6, Section 5] for a more recent survey with further references. For this purpose, we shall index each cluster by a finite sequence of positive integers, that is, by a vertex u of the Ulam–Harris tree ${\mathcal U}=\bigcup_{n=0}^{\infty} \mathbb{N}^n$ , such that the length $|u|$ of u corresponds to the generation of the cluster. By convention, the empty sequence $\varnothing$ with length 0 is used to label the cluster containing the ancestor. Clusters at the first generation are those such that there is a single untraceable edge along the contamination path from an infected individual in this cluster to the ancestor. They are indexed by $\mathbb{N}^1= \mathbb{N}=\{1,2, \ldots\}$ according to increasing order of their birth times, that is, times at which an individual in the ancestral cluster causes an untraceable contamination and generates a new cluster. For the sake of definitiveness, we agree that when the ancestral cluster generates only k untraceable contaminations until it is isolated, then the clusters indexed by $k+1, k+2, \ldots$ are fictitious clusters born at time $\infty$ . This is only a formality, and of course we shall only be concerned with non-fictitious clusters. By an obvious iteration, we label clusters at the nth generation by $u=(u_1, \ldots, u_n)\in \mathbb{N}^n$ for any $n\geq 0$ . It should be plain that the genealogy of clusters does not change with time, in the sense that once a cluster is born, its label will remain the same in the future.

For any $u\in {\mathcal U}$ , if the cluster labeled by u is not fictitious, then we write $\zeta_u$ for the age (that is, the time elapsed from the birth time) at which this cluster is isolated. We also write $\xi_u$ for the simple point process on $[0,\infty)$ of the ages at which this cluster is involved with untraceable contaminations, that is, it generates new clusters. So $\xi_u([0,t])$ is the number of children clusters generated when the cluster reaches age t, and in particular, $\xi_u([\zeta_u,\infty))=0$ . Finally, we write $C_u=(C_u(t), t\geq 0)$ , where $C_u(t)$ is the size of this cluster at time $t<\zeta_u$ , and we agree that $C_u(t)=0$ whenever $t\geq \zeta_u$ . In words, $C_u$ is the process of the number of contagious individuals in that cluster as a function of its age (recall that infected individuals are no longer contagious once they have been isolated); in particular, $C(\zeta\!-\!)$ is the size reached by the cluster at the time when it is detected. The most relevant information about the evolution of clusters and hence about the epidemic is encoded by the family of pairs

\begin{align*}{\mathbf C}_u=(C_u, \xi_u), \quad u\in \mathcal U,\end{align*}

where we agree for definitiveness that $C_u\equiv 0$ and $\xi_u\equiv 0$ when the cluster indexed by u is fictitious. Of course, ${\mathbf C}_u$ does not enable us to recover the subtree structure of the cluster indexed by u, but this is irrelevant for the questions we are interested in.

It should be intuitively clear that the distribution of the ancestral cluster ${\mathbf C}_{\varnothing}$ determines that of the whole process $({\mathbf C}_u)_{u\in \mathcal U}$ . More precisely, let us write ${\mathbf C}=(C, \xi)$ for a pair distributed as ${\mathbf C}_{\varnothing}$ , which we think of as describing the evolution of a typical cluster. Then it is readily checked that conditionally on $\xi_{\varnothing}(\mathbb{R}_+)=k$ , ${\mathbf C}_1, \ldots, {\mathbf C}_{k}$ are k independent and identically distributed (i.i.d.) copies of ${\mathbf C}$ that are further independent of ${\mathbf C}_{\varnothing}$ . More generally, it follows by iteration that for every $n\geq 1$ and any $u^1, \ldots, u^k\in \mathbb{N}^{n}$ , conditionally on the event that none of the clusters ${\mathbf C}_{u^1}, \ldots, {\mathbf C}_{u^k}$ are fictitious (which is measurable with respect to the family $({\mathbf C}_v\colon |v|<n)$ ), ${\mathbf C}_{u^1}, \ldots, {\mathbf C}_{u^k}$ are k i.i.d. copies of ${\mathbf C}$ that are further independent of $({\mathbf C}_v\colon |v|<n)$ . In other words, $({\mathbf C}_u)_{u\in \mathcal U}$ generates a Crump–Mode–Jagers branching process where the evolution of typical elements is distributed as ${\mathbf C}$ .

3. Statistics of a typical cluster

We discuss here some basic statistics of the typical cluster ${\mathbf C}=(C, \xi)$ in terms of the parameters $(\gamma,p,\delta)$ of the model. Recall that the integer-valued process C is absorbed at 0 at the time

\begin{align*}\zeta=\inf\{t\geq 0\colon C(t)=0\}\end{align*}

when this cluster is detected and isolated, and that $\xi$ is the point process of times at which untraceable contaminations occur.

It is now convenient to set

(3.1) \begin{equation} \rho \,{:\!=}\, \delta +p\gamma,\end{equation}

and recall that a Yule process with rate $\rho>0$ refers to a pure birth process with birth rate $\rho k$ from any state $k\geq 1$ and started from 1.

Lemma 3.1. The process $(C(t), t\geq 0)$ has the same law as

\begin{align*}(\textbf{1}_{\left\{{Y(t)\leq G}\right\}} Y(t), t\geq 0),\end{align*}

where $Y=(Y(t), t\geq 0)$ is a Yule process with rate $\rho$ , and G is an independent geometric variable with success probability $\delta/\rho$ , i.e. with tail distribution function

\begin{align*}\mathbb{P}(G \gt k)= (1-\delta/\rho)^{k} , \quad k\geq 0.\end{align*}

Lemma 3.1 shows in particular that the size $C(\zeta\!-\!)$ of the typical isolated cluster has the geometric distribution with success probability $\delta/\rho$ . The one-dimensional marginal laws of the typical cluster size process as well as the joint distribution of the time of isolation $\zeta$ and $C(\zeta\!-\!)$ follow readily.

Corollary 3.1. For every $t\geq 0$ , we have

\begin{align*}\mathbb{P}(C(t)=k)= (1-\delta/\rho)^{k-1} (1-\textrm{e}^{-\rho t})^{k-1} \,\textrm{e}^{-\rho t}\quad\ \textit{for}\ k\geq 1\end{align*}

and

\begin{align*} \mathbb{P}(C(t)=0)= \mathbb{P}(\zeta \leq t )=1- \dfrac{\rho}{\rho +\delta(\textrm{e}^{\rho t}-1)}.\end{align*}

Furthermore, we also have

\begin{align*}\mathbb{P}(C(\zeta\!-\!)=k, t\geq \zeta)= \dfrac{ \delta}{\rho}(1-\delta/\rho)^{k-1} (1-\textrm{e}^{-\rho t})^k \quad\ \textit{for}\ k\geq 1.\end{align*}

Proof. It suffices to write for $k\geq 1$ that

\begin{align*}\mathbb{P}(C(t)=k) &= \mathbb{P}(Y(t)=k , G\geq k) \\[3pt] &= \mathbb{P}(Y(t)=k) \mathbb{P}(G \geq k),\end{align*}

and recall that Y(t) has the geometric distribution with success probability $\textrm{e}^{-\rho t}$ . Then summation for $k\geq 1$ yields the second formula of the statement. We get the third formula similarly, writing for $k\geq 1$ that

\begin{align*}\mathbb{P}(C(\zeta\!-\!)=k, t\geq \zeta)= \mathbb{P}(G=k, Y(t)>k)= \mathbb{P}(G = k) \mathbb{P}(Y(t)> k).\\[-36pt]\end{align*}

We next turn our attention to the point process $\xi$ at which new clusters are generated, and write

\begin{align*}Z_1\,{:\!=}\, \xi(\mathbb{R}_+)\end{align*}

for the total number of (non-fictitious) clusters that the typical cluster begets. Its distribution is obtained by a slight variation of the argument for Lemma 3.1, and this entails the criterion for extinction of the epidemic that was stated in the Introduction.

Lemma 3.2. The variable $1+Z_1$ follows the geometric distribution with success probability ${\delta}/{((1-p) \gamma + \delta)}$ . In particular, $Z_1\in L^r(\mathbb{P})$ for all $r\geq 1$ ,

\begin{align*} \mathbb{E}(Z_1)= (1-p) \gamma / \delta ,\end{align*}

and as a consequence, the total number of infected individuals is finite (in other words, the epidemic eventually ceases) almost surely if and only if

\begin{align*}(1-p)\gamma \leq \delta.\end{align*}

Last, we introduce the intensity measure $\mu$ of the point process $\xi$ :

\begin{align*}\mu(t) \,{:\!=}\, \mathbb{E}(\xi([0,t])), \quad t\geq 0.\end{align*}

Corollary 3.2. For every $t\geq 0$ , there is the identity

\begin{align*} \mu(t) = (1-p)\dfrac{\gamma}{\delta} \biggl(1- \dfrac{ 1}{1+\delta(\textrm{e}^{\rho t}-1)/\rho}\biggr).\end{align*}

Proof. Indeed, the conditional probability given the process C, that an untraceable contamination event occurs during the time interval $[t, t+\textrm{d} t]$ , equals $(1-p)\gamma C(t)\,\textrm{d} t$ , and as a consequence,

\begin{align*}\textrm{d} \mu(t) = (1-p)\gamma\mathbb{E}(C(t))\,\textrm{d} t. \end{align*}

We deduce from Corollary 3.1 that

\begin{align*}\mathbb{E}(C(t))= \dfrac{ \textrm{e}^{\rho t}}{(1+\delta(\textrm{e}^{\rho t}-1)/\rho)^2}.\end{align*}

The formula in the statement follows.

We note that letting $t\to \infty$ in Corollary 3.2 yields

\begin{align*}\mathbb{E}(\xi(\mathbb{R}_+))=(1-p)\gamma/\delta,\end{align*}

in agreement with Lemma 3.2.

4. The Malthusian behavior

We shall assume throughout this section that

(4.1) \begin{equation} \delta \lt (1-p)\gamma,\end{equation}

so that the epidemic survives with strictly positive probability. More precisely, one immediately deduces from Lemma 3.2 that the probability of extinction equals $\delta/((1-p)\gamma)$ , which is the smallest solution to the equation $\mathbb{E}(x^{Z_1})=x$ . We shall derive here the main results of this work, simply by specifying in our setting some fundamental results of Nerman [Reference Nerman16] on the asymptotic behavior of Crump–Mode–Jagers branching processes with random characteristics. We start by introducing some of the key actors in this framework.

Consider the Laplace transform of the intensity measure of untraceable contaminations for a typical cluster:

\begin{align*}{\mathcal L}(x)= \int_0^{\infty} \,\textrm{e}^{-xt} \textrm{d} \mu(t), \quad x\geq 0.\end{align*}

Since ${\mathcal L}(0)=(1-p)\gamma/\delta \gt1$ , the equation $ {\mathcal L}(x)=1$ possesses a unique solution $\alpha=\alpha(\gamma,p,\delta)\in(0,\infty)$ , called the Malthusian parameter. That is, thanks to Corollary 3.2,

\begin{align*}(1-p)\gamma \int_0^{\infty} \dfrac{ \textrm{e}^{(\rho-\alpha) t}}{(1+\delta(\textrm{e}^{\rho t}-1)/\rho)^2}\,\textrm{d} t = 1,\end{align*}

or equivalently, in a slightly simpler form, using the change of variables $x=\textrm{e}^{-\rho t}$ ,

(4.2) \begin{equation} (1-p)\gamma \rho \int_0^1 \dfrac{ x^{\alpha/\rho}}{((\rho-\delta) x + \delta)^2}\,\textrm{d} x = 1. \end{equation}

We further set

(4.3) \begin{equation} \beta = - {\mathcal L}'(\alpha)= (1-p)\gamma \int_0^{\infty} t \dfrac{ \textrm{e}^{(\rho-\alpha) t}}{(1+\delta(\textrm{e}^{\rho t}-1)/\rho)^2}\,\textrm{d} t;\end{equation}

plainly $\beta\in (0,\infty)$ .

Next, it is convenient to use the notation

\begin{align*}\langle m,f\rangle \,{:\!=}\, \sum_{n=1}^{\infty} f(n) m(n),\end{align*}

where $m=(m(n), n\in \mathbb{N})$ is a finite measure on $\mathbb{N}$ and $f\colon \mathbb{N}\to \mathbb{R}_+$ is a generic non-negative function. We introduce two important measures $m^a$ and $m^i$ , related to typical active and isolated clusters respectively, by

\begin{align*}\langle m^a,f\rangle =\int_0^{\infty} \,\textrm{e}^{-\alpha t} \mathbb{E}(f(C(t)), t<\zeta) \,\textrm{d} t \end{align*}

and

\begin{align*} \langle m^i,f\rangle = \int_0^{\infty} \,\textrm{e}^{-\alpha t} \mathbb{E}(f(C(\zeta\!-\!)), \zeta\leq t ) \,\textrm{d} t .\end{align*}

These two measures can be determined explicitly from Corollary 3.1, using the notation

\begin{align*}\textrm{B}(x,y)= \dfrac{\Gamma(x) \Gamma(y)}{\Gamma(x+y)}=\int_0^1 s^{x-1}(1-s)^{y-1}\textrm{d} s , \quad x,y>0,\end{align*}

for the beta function. Indeed, from the change of variables $\textrm{e}^{-\rho t} = s$ , we then obtain that, for every $k\geq 1$ ,

(4.4) \begin{align} m^a(k) &= (1-\delta/\rho)^{k-1} \int_0^{\infty} \,\textrm{e}^{-\alpha t} (1-\textrm{e}^{-\rho t})^{k-1} \,\textrm{e}^{-\rho t} \textrm{d} t \nonumber \\[3pt] &= \dfrac{1}{\rho} (1-\delta/\rho)^{k-1} {\textrm{B}} (1+\alpha/\rho, k)\end{align}

and

(4.5) \begin{align} m^i(k) &=\dfrac{\delta}{\rho} (1-\delta/\rho)^{k-1} \int_0^{\infty} \,\textrm{e}^{-\alpha t} (1-\textrm{e}^{-\rho t})^{k} \textrm{d} t\nonumber \\[3pt] &= \dfrac{\delta}{\rho^2} (1-\delta/\rho)^{k-1} {\textrm{B}} (\alpha/\rho, k+1). \end{align}

Finally, we introduce

\begin{align*}W_n=\sum_{u\in \mathbb{N}^n} \,\textrm{e}^{-\alpha \sigma_u}, \quad n\geq 0,\end{align*}

where $\sigma_u$ stands for the birth time of the cluster labeled by u (so that $\sigma_u=\infty$ and $\textrm{e}^{-\alpha \sigma_u}=0$ if this cluster is fictitious). The process $(W_n, n\geq 0)$ is a martingale, often referred to as the intrinsic martingale; see Jagers [Reference Jagers9, Chapter 6]). Using the inequality $W_1\leq \xi(\mathbb{R}_+)$ and Lemma 3.2, we see that

\begin{align*}\mathbb{E}(W_1^2)<\infty,\end{align*}

and the uniform integrability of the intrinsic martingale follows; see e.g. [Reference Jagers10, Theorem 6.1]. We furthermore recall that its terminal value $W_{\infty}$ is strictly positive on the event that the epidemic survives, and of course $W_{\infty}=0$ on the event that the epidemic eventually ceases.

For the sake of simplicity, we focus on a few natural statistics of the epidemic at time $t\geq 0$ . Given a generic non-negative function $f\colon \mathbb{N}\to \mathbb{R}_+$ , we agree implicitly that $f(0)=0$ and write

\begin{align*} A^f(t) = \sum f(C_u(t-\sigma_u)),\end{align*}

where the sum is taken over all vertices u in the Ulam–Harris tree $ \mathcal U$ , such that the cluster labeled by u is born at time $\sigma_u\leq t$ (note that only active clusters at time t contribute to the sum). Turning our attention to isolated clusters, we similarly write

\begin{align*} I^f(t) =\sum f(C_u(\zeta_u\!-\!))\textbf{1}_{\left\{{\sigma_u+\zeta_u\leq t}\right\}},\end{align*}

where the sum is taken over all clusters that are isolated at time t. In the Crump–Mode–Jagers terminology, $A^f$ and $I^f$ are known as the processes counted with the random characteristics

(4.6) \begin{equation}\phi^a\colon t\mapsto f(C(t)) \textbf{1}_{\left\{{t<\zeta}\right\}} \quad\text{and} \quad \phi^i\colon t\mapsto f(C(\zeta\!-\!)) \textbf{1}_{\left\{{t\geq\zeta}\right\}},\end{equation}

respectively.

Recall the notation above, and notably (4.2), (4.3), (4.4), and (4.5). We can now state the main result of this work.

Theorem 4.1. Assume (4.1) and let $f\colon \mathbb{N}\to \mathbb{R}_+$ with $f(n)=O(\textrm{e}^{b n})$ for some $b\lt -\log(1-\delta/\rho)$ . The following limits then hold almost surely and in $L^1(\mathbb{P})$ :

\begin{align*}\lim_{t\to \infty} \,\textrm{e}^{-\alpha t} A^f(t) = \beta^{-1} \langle m^a, f \rangle W_{\infty}\end{align*}

and

\begin{align*} \lim_{t\to \infty} \,\textrm{e}^{-\alpha t} I^f(t)= \beta^{-1} \langle m^i, f \rangle W_{\infty}.\end{align*}

In particular, taking $f(n)=n$ yields the first-order asymptotic behavior of the total number of contagious (respectively isolated) individuals as time goes to infinity.

Proof. The claim of almost sure convergence is seen from Theorem 5.4 in Nerman [Reference Nerman16]; we just need to verify Conditions 5.1 and 5.2 there. For the first, we simply write

\begin{align*}\int_t^{\infty} \,\textrm{e}^{-\alpha s} \xi(\textrm{d} s) \leq \textrm{e}^{-\alpha t} \xi(\mathbb{R}_+),\end{align*}

and recall from Lemma 3.2 that $Z_1=\xi(\mathbb{R}_+)$ is integrable. This ensures that

\begin{align*}\mathbb{E}\biggl( \sup_{t\geq 0} \,\textrm{e}^{\alpha t} \int_t^{\infty} \,\textrm{e}^{-\alpha s} \xi(\textrm{d} s) \biggr) \lt \infty.\end{align*}

For the second, we assume for simplicity that $|f(n) |\leq \textrm{e}^{b n}$ for all $n\geq 1$ without loss of generality. The random characteristics in (4.6) can be bounded for all $t\geq 0$ by

\begin{align*}|\phi^a(t)| \leq \exp(b C(\zeta\!-\!)) \quad \text{and} \quad |\phi^i(t)| \leq \exp(b C(\zeta\!-\!)).\end{align*}

Observe that

\begin{align*}\mathbb{E}( \exp(b C(\zeta\!-\!))) \lt \infty,\end{align*}

since we know from Lemma 3.1 that $C(\zeta\!-\!)$ has the geometric distribution with success probability $\delta/\rho$ , and $(1-\delta/\rho) \,\textrm{e}^b \lt 1$ . It follows immediately that Condition 5.2 of [Reference Nerman16] holds for both $\phi^a$ and $\phi^i$ .

We next turn our attention to convergence in $L^1(\mathbb{P})$ . This in turn relies on [Reference Nerman16, Corollary 3.3], and we have to check equations 3.1 and 3.2 therein. The latter are both straightforward from the bounds established in the first part of this proof.

We next turn our attention to the empirical distribution of the sizes of active, respectively isolated, clusters. We denote the function identical to 1 on $\mathbb{N}$ by $\mathbf 1$ , so that in the above notation, $A^{\mathbf 1}(t)$ and $I^{\mathbf 1}(t)$ , respectively, are the number of active and of isolated clusters at time t. We then define the empirical distributions $\Pi^a(t)$ and $\Pi^i(t)$ for a generic function $f\colon \mathbb{N}\to \mathbb{R}_+$ by

\begin{align*}\langle \Pi^a(t), f\rangle = A^f(t)/A^{\mathbf 1}(t)\end{align*}

and

\begin{align*}\langle \Pi^i(t), f\rangle = I^f(t)/I^{\mathbf 1}(t).\end{align*}

We also introduce the normalized probability measures on $\mathbb{N}$ :

\begin{align*}\pi^a = m^a/ \langle m^a, \mathbf 1\rangle \quad \text{and} \quad \pi^i= m^i / \langle m^i, \mathbf 1 \rangle.\end{align*}

Thanks to (4.4) and (4.5), these are given explicitly by

(4.7) \begin{equation} \pi^a(k) = c_a(1-\delta/\rho)^{k-1} {\textrm{B}} (1+\alpha/\rho, k) \quad \text{for all }k\geq 1,\end{equation}

with

\begin{align*}1/c_a= \sum_{j=1}^{\infty} (1-\delta/\rho)^{j-1} {\textrm{B}} (1+\alpha/\rho, j),\end{align*}

and

(4.8) \begin{equation} \pi^i(k) = c_i(1-\delta/\rho)^{k-1} {\textrm{B}} (\alpha/\rho, k+1) \quad \text{for all }k\geq 1, \end{equation}

with

\begin{align*}1/c_i=\sum_{j=1}^{\infty} (1-\delta/\rho)^{j-1} {\textrm{B}} (\alpha/\rho, j+1).\end{align*}

We can now state the convergence of the empirical distributions.

Corollary 4.1. Assume (4.1). Then, conditionally on the event that the epidemic survives forever, we have almost surely

\begin{align*}\lim_{t\to \infty} \Pi^a(t) = \pi^a \quad\text{and} \quad \lim_{t\to \infty} \Pi^i(t) = \pi^i.\end{align*}

Proof. Indeed, recall that $W_{\infty}>0$ a.s. conditionally on survival of the epidemic. We derive from Theorem 4.1 that, on this event,

\begin{align*}\lim_{t\to \infty} A^f(t)/A^{\mathbf 1}(t) = \langle \pi^a, f\rangle\quad\text{and} \quad \lim_{t\to \infty} I^f(t)/I^{\mathbf 1}(t) = \langle \pi^i, f\rangle\end{align*}

for every bounded function $f\colon \mathbb{N}\to \mathbb{R}$ .

In words, Corollary 4.1 states that conditionally on survival of the epidemic, the empirical distributions of active cluster sizes and of isolated clusters sizes converge to $\pi^a$ and $\pi^i$ , respectively, as time goes to infinity. We shall therefore think of the latter as describing asymptotically the average distributions of the sizes of active clusters and isolated clusters, respectively. It is interesting to point out similarities between Corollary 4.1 and earlier results by Deijfen [Reference Deijfen5, Theorem 1.1 and Example 2] on the asymptotic degree distribution for certain random evolving networks. More specifically, in this model new vertices arrive in continuous time, are connected to an existing vertex with probability proportional to the so-called fitness of that vertex, and vertices then die at rates depending on their accumulated in-degrees. Although the model considered by Deijfen is different from ours, it also bears a clear resemblance, and the similarities between the results (and the methods as well) should not come as a surprise.

It is also interesting to observe from formulas (4.7) and (4.8) and the elementary identity

\begin{align*}\dfrac{\alpha}{\rho} {\textrm{B}} (\alpha/\rho, k+1) = k {\textrm{B}} (1+\alpha/\rho, k)\end{align*}

that

\begin{align*}\pi^i(k) = \dfrac{k\pi^a(k)}{\sum_{j=1}^{\infty} j \pi^a(j)} \quad \text{for all }k\geq 1.\end{align*}

In words, the average distribution of the sizes of isolated clusters is the size-biased version of that of active clusters. This relation stems from the fact that the rate at which an active cluster becomes isolated is proportional to its size. Since the empirical distribution of active cluster sizes converges to $\pi^a$ , the empirical distribution of isolated cluster sizes must converge to the size-biased version of $\pi^a$ . We refer to Corollary 1 of [Reference Bansaye, Gu and Yuan2] and its proof for details of a rigorous argument.

5. Concluding comments

5.1. Comparison with a model of Bansaye, Gu, and Yuan, and an eigenproblem

This work was inspired by a recent manuscript of Bansaye et al. [Reference Bansaye, Gu and Yuan2], in which they introduced a similar model for epidemics with contact tracing and cluster isolation. The main difference from the present article is that in [Reference Bansaye, Gu and Yuan2], contaminations are always traceable initially, but traceability gets lost at some fixed rate. In other words, edges in the contamination tree are traceable when they first appear, and become untraceable after exponentially distributed time-laps, independently of the other edges. As a consequence, clusters do not only grow when new contamination events occur, but also split when a traceable edge becomes untraceable.

Bansaye et al. investigate the large time asymptotic behavior of the epidemic using different tools, namely they first analyze a deterministic eigenproblem for a growth–fragmentation–isolation equation that is naturally related to their setting; furthermore they also rely on known properties of random recursive trees. They establish results similar to our Theorem 4.1 and Corollary 4.1 in terms of these eigenelements; the statements in [Reference Bansaye, Gu and Yuan2] are, however, less precise than ours, as no explicit formulas for the eigenelements are given (only their existence is established).

In our setting, using the notation of Section 4, the expectation of linear functionals of clusters at a given time yields a family $(\nu_t, t\geq 0)$ of measures on $\mathbb{N}$ given by

\begin{align*}\langle \nu_t, f\rangle = \mathbb{E}(A^f(t)),\end{align*}

where $f\colon \mathbb{N}\to \mathbb{R}_+$ is a generic bounded function. From the dynamics of the epidemic, we get the evolution equation

(5.1) \begin{equation} \textrm{d} \langle \nu_t, f\rangle = \langle \nu_t, {\mathcal A}f\rangle \,\textrm{d} t, \end{equation}

with

(5.2) \begin{equation} {\mathcal A}f(k) =k( p\gamma (f(k+1)-f(k)) + (1-p)\gamma f(1)- \delta f(k)); \end{equation}

the initial condition $\nu_0$ is the Dirac mass at 1 since we assume that the epidemic starts from a single contagious individual. Specifically, in (5.2), the term $kp\gamma (f(k+1)-f(k))$ accounts for the growth of a cluster from size k to $k+1$ , which occurs with rate $kp\gamma$ . The term $k(1-p)\gamma f(1)$ stems from the birth of new clusters of size 1 (i.e. an untraceable contamination) induced by a cluster of size k, which occurs with rate $k(1-p)\gamma$ , and finally $-k\delta f(k)$ for the isolation of a cluster of size k, which occurs with rate $k\delta$ . This formula for the infinitesimal generator $ {\mathcal A}$ should be compared with Lemma 1 in [Reference Bansaye, Gu and Yuan2], and notably equation (4.15) therein.

Predominantly, growth–fragmentation equations (and more generally, evolution equations) cannot be solved explicitly, and most works in this area are concerned with the large time asymptotic behavior of its solutions; see [Reference Bansaye, Gu and Yuan2] for some references. Roughly speaking, the paradigm, which stems from the Perron–Frobenius theorem for matrices with positive entries, is to resolve the eigenproblem for the infinitesimal generator, that is, to determine the principal eigenvalue (i.e. the eigenvalue with the largest real part) and its left eigenfunctions. The principal eigenvalue is identified as the Malthusian parameter, and the left eigenfunction (viewed as a measure) yields the so-called asymptotic profile, that is, in our setting, the measure $m^a$ in Theorem 4.1. So the analysis carried out in the present Section 4 solves this eigenproblem indirectly for (5.2), the solution being given by (4.2) and (4.7). Specifically, it holds for all bounded $f\colon \mathbb{N}\to \mathbb{R}_+$ that

\begin{align*}\langle \mathcal{A}^{\top}\pi^a , f \rangle \,{:\!=}\,\langle \pi^a ,\mathcal{A} f \rangle = \alpha\langle \pi^a, f\rangle.\end{align*}

However, it does not seem straightforward to check this identity directly, and we shall provide more details below.

Let $\nu_t(k)$ denote the expected number of active clusters of size $k\geq 1$ at time t. From the dynamics of the epidemic (see the discussion following (5.2)), we have that for $k\geq 2$

(5.3) \begin{equation} \dfrac{\partial \nu_t(k)}{\partial t} + p\gamma( k \nu_t(k)- (k-1)\nu_t(k-1)) =- \delta k \nu_t(k) ,\end{equation}

whereas for $k=1$

(5.4) \begin{equation} \dfrac{\partial \nu_t(1)}{\partial t} + p\gamma \nu_t(1) = (1-p)\gamma \sum_{j=1}^{\infty} j\nu_t(j) -\delta \nu_t(1). \end{equation}

Of course, (5.3) and (5.4) are equivalent to the evolution equation (5.1). From the point of view of age-structured population models (recall Remark 2.1), these should be viewed as a version of the McKendrick–von Foerster PDE; see [Reference Kot12, equations (23.4) and (23.5)].

Following [Reference Kot12, Chapter 23], it is then natural to search for a special solution to (5.3) and (5.4) in the form $\nu_t(k)= \textrm{e}^{r t}\nu(k)$ for some $r>0$ and some measure $\nu$ on $\mathbb{N}$ , which of course amounts to solving the eigenproblem $ \mathcal A^{\top} \nu = r \nu$ . Recall the notation (3.1); from (5.3) and (5.4), we first get the linear recurrence equation

(5.5) \begin{equation} \nu(k) = \dfrac{p\gamma (k-1)}{r+\rho k}\nu(k-1), \quad k\geq 2, \end{equation}

and then, for $k=1$ , the identity

(5.6) \begin{equation} (r+\rho) \nu(1) = (1-p)\gamma \sum_{j=1}^{\infty} j\nu(j).\end{equation}

We readily deduce from (5.5) and well-known properties of the beta function $\textrm{B}$ that

(5.7) \begin{equation} \nu(k) = c(1-\delta/\rho)^{k-1} \textrm{B}(1+r/\rho,k ), \quad k\geq 1,\end{equation}

where $c>0$ is some arbitrary constant. We note that $(r+\rho)\nu(1)/c=\rho$ , and can now determine r by rewriting (5.6) in the form

\begin{align*}\rho &= (1-p)\gamma \sum_{j=1}^{\infty} j (1-\delta/\rho)^{j-1} \textrm{B}(1+r/\rho,j )\\[3pt] & = (1-p)\gamma \sum_{j=1}^{\infty} \int_0^1 j (1-\delta/\rho)^{j-1} (1-x)^{j-1} x^{r/\rho}\textrm{d} x \\[3pt] & = (1-p)\gamma \int_0^1 \dfrac{ x^{r/\rho}}{(1-(1-\delta/\rho)(1-x))^2}\,\textrm{d} x \\[3pt] & = (1-p)\gamma \rho^2 \int_0^1 \dfrac{ x^{r/\rho}}{((\rho-\delta)x+\delta)^2}\,\textrm{d} x. \end{align*}

We have recovered (4.2), which determines the Malthusian parameter. So $r=\alpha$ , and we conclude from (5.7) and (4.7) that $\nu$ and $\pi^a$ are indeed proportional.

The calculations above are reminiscent of those for the Leslie model [Reference Kot12, Chapter 22], and in particular (4.2) can be thought of as an Euler–Lotka equation [Reference Kot12, equation (20.6)].

5.2. A detection paradox

We next discuss in more detail the detection paradox mentioned in the Introduction. Imagine that we rank the clusters in increasing order of their birth times rather than indexing them by the Ulam–Harris genealogical tree as we did previously. This sequence is infinite if and only if the epidemic survives, and conditionally on that event, its elements are independent, each being distributed as the typical cluster. One may then expect from the law of large numbers that as time goes to infinity, the limit $\pi^i$ of the empirical distribution of the sizes of isolated clusters should coincide with the law of the size of a typical cluster at the time when it is detected, that is, by Lemma 3.1, the geometric distribution with success probability $\delta/\rho$ . However, (4.8) shows that this is not the case, and more precisely, $\pi^i$ is a biased version of the geometric law, where the bias is given by a beta function.

The naive argument above of course has a flaw, which stems from the fact that the empirical distribution of the isolated clusters at a given time corresponds to a partial sum of clusters which are listed in increasing order of their detection times rather than their birth times. This reordering tends to list first clusters which are quickly detected and hence had little time to grow, which hints at the feature that the average isolated cluster size is dominated stochastically by the size of a typical isolated cluster. Nonetheless, reordering alone is not sufficient to explain the detection paradox: the second crucial ingredient is the exponential growth, and more precisely the fact that the number of clusters, say for simplicity born during the time interval $[t,t+1]$ , is of the same order as the number of all the clusters born before time t, no matter how large t is. A significant proportion of clusters born during $[t,t+1]$ are detected before time $t+1$ ; due to the time constraint, these clusters have on average a smaller size than the typical cluster when it is isolated, and this explains the seeming paradox.

For a better understanding of the mechanisms at work in the explanation above, it may be useful to consider the following elementary example. Consider a Poisson point process on $\mathbb{R}_+\times (0,\infty)$ whose atoms are denoted generically by $(b,\ell)$ , and which has intensity $\textrm{e}^{b} \textrm{d} b \lambda(\textrm{d} \ell)$ , where $\lambda$ is some probability measure on $(0,\infty)$ . We think of $(b, \ell)$ as an individual born at time b and with lifespan $\ell$ . Imagine that we want to estimate the lifespan distribution $\lambda$ , that is, more specifically, the quantity $\langle \lambda,f\rangle$ for an arbitrary bounded continuous function $f\colon (0,\infty)\to \mathbb{R}$ , from the observation of the population up to some large time t. If we could observe the lifespan of an individual at the time when it is born, then this would be an easy matter. Indeed, it then suffices to compute the empirical mean of $f(\ell)$ for individuals $(b,\ell)$ born at time $b\leq t$ , and it is readily checked by Poissonian computation that this quantity converges almost surely to $\langle \lambda,f\rangle$ as $t\to \infty$ . But of course it is unrealistic to assume that the lifespan can be observed at the birth of an individual, and let us instead assume that lifespan can be observed at death only.

The total number of dead individuals at time t has the Poisson distribution with parameter

\begin{align*}\int_0^t \,\textrm{e}^b \lambda((0,t-b]) \,\textrm{d} b \sim \textrm{e}^t \int_{(0,\infty)} \,\textrm{e}^{-\ell} \lambda(\textrm{d} \ell) \quad \text{as $t\to \infty$.} \end{align*}

More generally, it is easily checked that if we write $\langle M(t),f\rangle$ for the empirical mean of $f(\ell)$ computed for all individuals who are dead at time t, that is, such that $b+\ell\leq t$ , then

\begin{align*}\lim_{t\to \infty} \langle M(t),f\rangle =\langle \lambda_1,f\rangle,\end{align*}

where

\begin{align*}\lambda_1(\textrm{d} \ell) = \dfrac{\textrm{e}^{-\ell} \lambda(\textrm{d} \ell)}{\int_{(0,\infty)} \,\textrm{e}^{-s} \lambda(\textrm{d} s)}.\end{align*}

In other words, the empirical mean $\langle M(t),f\rangle$ is a consistent estimator of $ \langle \lambda_1,f\rangle$ rather than of $\langle \lambda,f\rangle$ .

We stress that this detection paradox disappears for a version of this model where the intensity of the Poisson point process only grows sub-exponentially in time – say for simplicity it is given by $b^r \textrm{d} b \lambda(\textrm{d} \ell)$ for some $r>0$ . The same calculation as above easily shows that the empirical mean of $f(\ell)$ computed for all individuals who are dead by time t then does converge to $\langle \lambda,f\rangle$ as $t\to \infty$ .

5.3. Relation to the Yule–Simon distribution

In 1955, following G. Udny Yule [Reference Yule19], Herbert A. Simon [Reference Simon18] introduced an elementary model depending on a parameter $q\in(0,1)$ (that accounts for the memory of the model), which today would be referred to as an algorithm with preferential attachment. Simon’s algorithm produces a random text, that is, a long string of words $w_1\ldots w_n$ , as follows. Once the first word $w_1$ has been written, for each $j=1, \ldots,n-1$ , $w_{j+1}$ is copied from a uniform sample from $w_1, \ldots, w_j$ with probability q, and with complementary probability $1-q$ , $w_{j+1}$ is a new word different from all the preceding words. Simon proved that for every fixed $k\geq 1$ , the expected proportion of different words that have been written exactly k times in the text converges as $n\to \infty$ towards

(5.8) \begin{equation} \sigma_q(k)=\dfrac{1}{q} \textrm{B}(1+1/q, k).\end{equation}

The probability measure on $\mathbb{N}$ , $\sigma_q=(\sigma_q(k), k\geq 1)$ , is known as the Yule–Simon distribution with parameter $1/q$ . Comparing (4.7) with (5.8), we can now view the average distributions of the sizes of active clusters $m^a$ as an exponentially tilted version of the Yule–Simon distribution with parameter $\alpha/\rho$ .

In this direction, we observe that the limiting case of our model with $\delta=0$ , which corresponds to a degenerate case where detection is absent, merely rephrases Simon’s algorithm with memory parameter $q=p$ . When there is no detection, the evolution of a typical cluster is just that of a Yule process with rate $p\gamma$ (without killing). Then the intensity measure of birth times of new clusters given by $\mu(\textrm{d} t)= (1-p) \gamma \,\textrm{e}^{p\gamma t} \textrm{d} t$ and the Malthusian parameter can be identified by solving

\begin{align*}(1-p) \gamma \int_0^{\infty} \,\textrm{e}^{(p\gamma-\alpha)t}\textrm{d} t=1.\end{align*}

Plainly we have $\alpha = \gamma$ and the parameter of the Yule–Simon distribution is simply $\alpha/\rho=1/p$ . In this setting, the degenerate case of Corollary 4.1 for $\delta =0$ can be viewed as a strong version of Simon’s result, where the convergence is almost sure and not just for the expectation. See [Reference Holmgren and Janson6, Example B.11] for a closely related discussion in the setting of Yule’s original model of evolution of species, which is a bit different but nonetheless also yields the Yule–Simon distribution (5.8).