Branching processes in random environments with thresholds

1. Introduction

Branching processes and their variants are used to model various biological, biochemical, and epidemic processes [Reference Jagers1–Reference Kimmel and Axelrod4]. More recently, these methods have been used to model the spread of COVID cases in communities during the early stages of the pandemic [Reference Yanev, Stoimenova and Atanasov5, Reference Atanasov, Stoimenova and Yanev6]. As time progressed, varying local containment efforts caused changes in the number of infected members in each community [Reference Falcó and Corral7, Reference Sun, Kryven and Bianconi8], leading to periods of increase and decrease. In this paper, we describe a stochastic process model built on a branching process in random environments (BPRE) that explicitly takes into account periods of growth and decrease in the transmission rate of the virus.

Specifically, we consider a branching process model initiated by a random number of ancestors (thought of as initiators of the pandemic within a community). During the first several generations, the process grows uncontrolled, allowing immigration into the system. This initial phase is modeled using a supercritical branching process with immigration in random environments, specifically independent and identically distributed (i.i.d.) environments. When consequences of rapid spread become significant, policymakers introduce restrictions to reduce the rate of growth, hopefully resulting in a reduced number of infected cases. The limitations are modeled using upper thresholds on the number of infected cases, and beyond the threshold the process changes its character to evolve as a subcritical branching process in random environments. During this period—owing to strict controls—immigration is also not allowed. In practical terms, this period typically involves a ‘lockdown’ and other social containment efforts, the intensity of which varies across communities.

The period of restrictions is not sustainable for various reasons, including political, social, and economic pressures leading to the easing of controls. Policymakers use multiple metrics to gradually reduce controls, leading to an ‘opening of communities’, resulting in increased human interaction. As a result, or because of changes undergone by the virus, the number of infected cases increases again. We use lower thresholds in the number of ‘newly infected’ to model the period of change and let the process evolve again as a supercritical BPRE in i.i.d. environments after it crosses the lower threshold. The process continues to evolve in this manner, alternating between periods of increase and decrease. In this paper, we provide a rigorous probabilistic analysis of this model.

Although we have taken the dynamics of COVID spread as a motivation for the proposed model, the aforementioned cyclic behavior is often observed in other biological systems, such as those modeled by predator–prey models or the susceptible–infected–recovered (SIR) model. In some biological populations, the cyclical behavior can be attributed to a decline in fecundity as the population size approaches some threshold [Reference Klebaner9]. Deterministic models such as ordinary differential equations, dynamical systems, and corresponding discrete-time models are used for analysis in the applications mentioned above [Reference Teschl10–Reference Iannelli and Pugliese12]. While many of the models described above yield good qualitative descriptions, uncertainty estimates are typically unavailable. It is worth pointing out that the previously described branching process methods also produce reasonable point estimates for the mean growth during the early stages of the pandemic. However, these point estimates are unreliable during the later stages of the pandemic. In this paper, we address statistical estimation of the mean growth and characterize the variance of the estimates. We end the discussion with a plot, Figure 1, of the total number of confirmed COVID cases per week in Italy from 23 February 2020 to 20 July 2022. The plot also includes the number of cases estimated using the proposed model. Other examples with similar plots include the hare–lynx predator–prey dynamics and measles cases [Reference Iannelli and Pugliese12–Reference Hempel and Earn14].

Figure 1. In black weekly COVID cases in Italy from February 23, 2020 to February 3, 2023. In blue a BPRE starting with the same initial value and offspring mean having the negative binomial distribution with predefined number of successful trials $r=10$ and Gamma-distributed mean with shape parameter equal to the mean of the data and rate parameter 1.

Before we provide a precise description of our model, we begin with a brief description of BPREs with immigration. Let $\Pi_{n} = \big(P_{n},Q_{n}\big)$ be i.i.d. random variables taking values in $\mathcal{P} \times \mathcal{P}$ , where $\mathcal{P}$ is the space of probability distributions on $\mathbb{N}_{0}$ ; that is, $P_{n}=\{ P_{n,r} \}_{r=0}^{\infty}$ and $Q_{n}=\{ Q_{n,r} \}_{r=0}^{\infty}$ for some non-negative integers $P_{n,r}$ and $Q_{n,r}$ such that $\sum_{r=0}^{\infty} P_{n,r}=1$ and $\sum_{r=0}^{\infty} Q_{n,r}=1$ . The process $\Pi=\{ \Pi_{n} \}_{n=0}^{\infty}$ is referred to as the environmental sequence. For each realization of $\Pi$ , we associate a population process $\{ Z_{n} \}_{n=0}^{\infty}$ defined recursively as follows: let $Z_{0}$ take values on the positive integers, and for $n \geq 0$ , let

\begin{equation*} Z_{n+1}=\sum_{i=1}^{Z_{n}} \xi_{n,i}+I_{n},\end{equation*}

where, given $\Pi_{n}=\big(P_{n},Q_{n}\big)$ , $\{ \xi_{n,i} \}_{i=1}^{\infty}$ are i.i.d. with distribution $P_{n}$ and $I_{n}$ is an independent random variable with distribution $Q_{n}$ . The random variable $Y_{n}=\log(\overline{P}_{n})$ , where $\overline{P}_{n}=\sum_{r=0}^{\infty} r P_{n,r}$ , plays an important role in the classification of BPREs with immigration. It is well known that when $\mathbb{E}[Y_{0}] > 0$ , the process diverges to infinity with probability one, and if $\mathbb{E}[Y_{0}] \leq 0$ and the immigration is degenerate at zero for all environments, then the process becomes extinct with probability one [Reference Athreya and Karlin15]. Furthermore, in the subcritical case, that is, $\mathbb{E}[Y_{0}] < 0$ , one can identify three distinct regimes: (i) weakly subcritical, (ii) moderately subcritical, and (iii) strongly subcritical. The regime (i) corresponds to the case when there exists a $0 < \rho < 1$ such that $\mathbb{E}\big[Y_{0} e^{\rho Y_{0}}\big]=0$ , while (ii) corresponds to the case when $\mathbb{E}[Y_{0} e^{Y_{0}}]=0$ . Finally, (iii) corresponds to the case when $\mathbb{E}[Y_{0} e^{Y_{0}}]<0$ [Reference Kersting and Vatutin16]. In this paper, when working with the subcritical regime, we will assume that the process is strongly subcritical, and we will refer to it as a subcritical process in the rest of the manuscript.

We now turn to a description of the model. Let $\Pi^{U} = \{ \Pi_{n}^{U} \}_{n=0}^{\infty}$ , where $\Pi_{n}^{U}=\big(P_{n}^{U},Q_{n}^{U}\big)$ , denote a collection of supercritical environmental sequences. Here, $P_{n}^{U}=\{ P_{n,r}^{U} \}_{r=0}^{\infty}$ indicates the offspring distribution and $Q_{n}^{U}=\{ Q_{n,r}^{U} \}_{r=0}^{\infty}$ represents the immigration distribution. Also, let $\Pi^{L}=\{\Pi_{n}^{L}\}_{n=0}^{\infty}$ , where $\Pi_{n}^{L}=P_{n}^{L}=\{ P_{n,r}^{L} \}_{r=0}^{\infty}$ , denote a collection of subcritical environmental sequences. We now provide an evolutionary description of the process: at time zero the process starts with a random number of ancestors $Z_{0}$ . Each of them lives one unit of time and reproduces according to the distribution $P_{0}^{U}$ . Thus, the size of the first-generation population is

\begin{equation*} Z_{1}=\sum_{i=1}^{Z_{0}} \xi_{0,i}^{U}+I_{0}^{U},\end{equation*}

where, given $\Pi_{0}^{U}=(P_{0}^{U},Q_{0}^{U})$ , the $\xi_{0,i}^{U}$ are i.i.d. random variables with offspring distribution $P_{0}^{U}$ and are independent of the immigration random variable $I_{0}^{U}$ with distribution $Q_{0}^{U}$ . The random variable $\xi_{0,i}^{U}$ is interpreted as the number of children produced by the ith parent in the 0th generation, and $I_{0}^{U}$ is interpreted as the number of immigrants whose distribution is generated by the same environmental random variable $\Pi^{U}_{0}$ .

Let $U_{1}$ denote the random variable representing the upper threshold. If $Z_{1} < U_{1}$ , each member of the first-generation population lives one unit of time and evolves, conditionally on the environment, as the ancestors independent of the population size at time one. That is,

\begin{equation*} Z_{2} = \sum_{i=1}^{Z_{1}} \xi_{1,i}^{U}+I_{1}^{U}.\end{equation*}

As before, given $\Pi_{1}^{U}=\big(P_{1}^{U},Q_{1}^{U}\big)$ , the $\xi_{1,i}^{U}$ are i.i.d. with distribution $P_{1}^{U}$ and $I_{1}^{U}$ has distribution $Q_{1}^{U}$ . The random variables $\xi_{1,i}^{U}$ are independent of $Z_{1}$ , $\xi_{0,i}^{U}$ , and $I_{0}^{U}$ , $I_{1}^{U}$ . If $Z_{1} \geq U_{1}$ , then

\begin{equation*} Z_{2} = \sum_{i=1}^{Z_{1}} \xi_{1,i}^{L},\end{equation*}

where, given $\Pi_{1}^{L}=P_{1}^{L}$ , the $\xi_{1,i}^{L}$ are i.i.d. with distribution $P_{1}^{L}$ . Thus, the size of the second-generation population is

\begin{equation*} Z_{2} = \begin{cases} \sum_{i=1}^{Z_{1}} \xi_{1,i}^{U}+I_{1}^{U} &\text{ if } Z_{1}<U_{1}, \\[5pt] \sum_{i=1}^{Z_{1}} \xi_{1,i}^{L} &\text{ if } Z_{1} \geq U_{1}. \end{cases}\end{equation*}

The process $Z_{3}$ is defined recursively as before. As an example, if $Z_{1} < U_{1}$ , $Z_{2} < U_{1}$ or $Z_{1} \geq U_{1}$ , $Z_{2} \leq L_{1}$ , for a random lower threshold $L_{1}$ , then the process will evolve like a supercritical BPRE with offspring distribution $P_{2}^{U}$ and immigration distribution $Q_{2}^{U}$ . Otherwise (that is, $Z_1< U_1$ and $Z_2 \ge U_1$ or $Z_1 \ge U_1$ and $Z_2 >L_1$ ), the process will evolve like a subcritical BPRE with offspring distribution $P_{2}^{L}$ . This dynamics continues with different thresholds $(U_j, L_j)$ , yielding the process $\{ Z_{n} \}_{n=0}^{\infty}$ , which we refer to as a branching process in random environments with thresholds (BPRET). The consecutive set of generations where the reproduction is governed by a supercritical BPRE is referred to as the supercritical regime, while the other is referred to as the subcritical regime. As we will see below, non-trivial immigration in the supercritical regime is required to obtain alternating periods of increase and decrease.

The model described above is related to size-dependent branching processes with a threshold as studied by Klebaner [Reference Klebaner9] and more recently by Athreya and Schuh [Reference Athreya and Schuh17]. Specifically, in that model the offspring distribution depends on a fixed threshold K and the size of the previous generation. As observed in these papers, these Markov processes either explode to infinity or are absorbed at zero. In our model the thresholds are random and dynamic, resulting in a non-Markov process; however, the offspring distribution does not depend on the size of the previous generation as long as they belong to the same regime. Indeed, when $U_{j}-1=L_{j}=K$ for all $j \geq 1$ , the immigration distribution is degenerate at zero, and the environment is fixed, one obtains as a special case the density-dependent branching process (see for example [Reference Klebaner9, Reference Athreya and Schuh17–Reference Jagers and Klebaner20]. Additionally, while the model of Klebaner [Reference Klebaner9] uses Galton–Watson processes as a building block, our model uses branching processes in i.i.d. environments.

Continuing with our discussion on the literature, Athreya and Schuh [Reference Athreya and Schuh17] show that in the fixed-environment case, the special case of a size-dependent process with a single threshold becomes extinct with probability one. We show that this is also the case for the BPRE when there is no immigration; the details are in Theorem 2.1. Similar phenomena have been observed in slightly different contexts in Jagers and Zuyev [Reference Jagers and Zuyev21, Reference Jagers and Zuyev22]. The incorporation of an immigration component ensures that the process is not absorbed at zero and hence may be useful for modeling stable populations at equilibrium as done in deterministic models. For additional discussion see Section 7.

For ease of further discussion, we introduce some notation. Let $Y_{n}^{U} \;:\!=\; \log\!\Big(\overline{P}_{n}^{U}\Big)$ and $Y_{n}^{L} \;:\!=\; \log\!\Big(\overline{P}_{n}^{L}\Big)$ , where

\begin{equation*} \overline{P}_{n}^{U}=\sum_{r=0}^{\infty} r P_{n,r}^{U} \quad \text{ and } \quad\overline{P}_{n}^{L}=\sum_{r=0}^{\infty} r P_{n,r}^{L};\end{equation*}

that is, $\overline{P}_{n}^{U}$ and $\overline{P}_{n}^{L}$ represent the offspring means conditional on the environments $\Pi_{n}^{U}=\big(P_{n}^{U},Q_{n}^{U}\big)$ and $\Pi_{n}^{L}=P_{n}^{L}$ , respectively. Also, let $\overline{Q}_{n}^{U}=\sum_{r=0}^{\infty} r Q_{n,r}^{U}$ denote the immigration mean conditional on the environment, and let

\begin{equation*} \overline{\overline{P}}_{n}^{U}=\sum_{r=0}^{\infty} \Big(r-\overline{P}_{n}^{U}\Big)^{2} P_{n,r}^{U} \quad \text{ and } \quad \overline{\overline{P}}_{n}^{L}=\sum_{r=0}^{\infty} \Big(r-\overline{P}_{n}^{L}\Big)^{2} P_{n,r}^{L}\end{equation*}

denote the conditional variance of the offspring distributions given the environment.

From the description, it is clear that the crossing times at the thresholds $(U_{j},L_{j})$ of $Z_{n}$ , namely $\tau_{j}$ and $\nu_{j}$ , will play a significant role in the analysis. It will turn out that $\{Z_{\tau_{j}}\}$ and $\{Z_{\nu_{j}}\}$ form time-homogeneous Markov chains with state spaces $S^{L} \;:\!=\; \mathbb{N}_{0} \cap [0,L_{U}]$ and $S^{U} \;:\!=\; \mathbb{N} \cap [L_{U}+1,\infty)$ , respectively, where we take $L_{j} \leq L_{U}$ and $U_{j} \geq L_{U}+1$ for all $j \geq 1$ . Under additional conditions on the offspring distribution and the environment sequence, the processes $\{Z_{\tau_{j}}\}$ and $\{Z_{\nu_{j}}\}$ will be uniformly ergodic. These results are established in Section 3.

The amount of time the process spends in the supercritical and subcritical regimes, beyond its mathematical and scientific interest, will also arise in the study of the central limit theorem for the estimates of $M^{U} \;:\!=\; \mathbb{E}\Big[\overline{P}_{n}^{U}\Big]$ and $M^{L} \;:\!=\; \mathbb{E}\Big[\overline{P}_{n}^{L}\Big]$ . Using the uniform ergodicity alluded to above, we will establish that the time averages of $\tau_{j}-\nu_{j-1}$ and $\nu_{j}-\tau_{j}$ converge to finite positive constants, $\mu^{U}$ and $\mu^{L}$ . Additionally, we establish a central limit theorem related to this convergence under a finite-second-moment hypothesis after an appropriate centering and scaling, that is,

\begin{align*} \frac{1}{\sqrt{n}}\sum_{j=1}^n (\tau_j -\nu_{j-1}) \xrightarrow[n \to \infty]{d} N\big( \mu^U, \sigma^{2,U}\big),\end{align*}

and we characterize $\sigma^{2,U}$ in terms of the stationary distribution of the Markov chain. A similar result also holds for $\nu_j-\tau_{j}$ . This, in turn, provides qualitative information regarding the proportion of time the process spends in these regimes. That is, if $C_{n}^{U}$ is the amount of time the process spends in the supercritical regime up to time $n-1$ , we show that $n^{-1}C_n^{U}$ converges to $\mu^U\big(\mu^U+\mu^L\big)^{-1}$ ; a related central limit theorem is also established, and in the process we characterize the limiting variance. Interestingly, we show that the central limit theorem prevails even for the joint distribution of the length of time and the proportion of time the process spends in the supercritical and subcritical regimes. These results are described in Sections 4 and 5.

An interesting question concerns the rate of growth of the BPRET in the supercritical and subcritical regimes described by the corresponding expectations, namely $M^{U}$ and $M^{L}$ . Specifically, we establish that the limiting joint distribution of the estimators is bivariate normal with a diagonal covariance matrix, yielding asymptotic independence of the mean estimators derived using data from supercritical and subcritical regimes. In the classical setting of a supercritical BPRE without immigration, this problem has received some attention (see for instance Dion and Esty [Reference Yanev, Stoimenova and Atanasov23]). The problem considered here is different in the following four ways: (i) the population size does not converge to infinity, (ii) the lengths of the regimes are random, (iii) in the supercritical regime the population size may be zero, and (iv) there is an additional immigration term. While (iii) and (iv) can be accounted for in the classical settings as well, their effect on the point estimates is minimized because of the exponential growth of the population size. Here, while the exponential growth is ruled out, perhaps as anticipated, the Markov property of the process at crossing times, namely $\{Z_{\tau_j}\}$ and $\{Z_{\nu_j}\}$ , and their associated regeneration times plays a central role in the proof. It is important to note that it is possible for both regimes to occur between regeneration times. Hence, the proportion of time that the process spends in the supercritical and subcritical regimes also plays a vital role in the derivation of the asymptotic limit distribution. The limiting variance of the estimators depends additionally on $\mu^{U}$ and $\mu^{L}$ , beyond $V_{1}^{U} \;:\!=\; \mathbb{V}\Big[\overline{P}_{0}^{U}\Big]$ , $V_{1}^{L} \;:\!=\; \mathbb{V}\Big[\overline{P}_{0}^{L}\Big]$ , $V_{2}^{U} \;:\!=\; \mathbb{E}\Big[\overline{\overline{P}}_{0}^{U}\Big]$ , and $V_{2}^{L} \;:\!=\; \mathbb{E}\Big[\overline{\overline{P}}_{0}^{L}\Big]$ . In the special case of fixed environments, the limit behavior of the estimators takes a different form compared to the traditional results, as described for example in Heyde [Reference Heyde24]. These results are in Section 6.

Finally, in Appendix B we provide some numerical experiments illustrating the behavior of the model. Specifically, we illustrate the effects of different distributions on the path behavior of the process and describe how they change when the thresholds increase. The experiments also suggests that if different regimes are not taken into account, the true growth rate of the virus may be underestimated. We now turn to Section 2, where we develop additional notation and provide precise statements of the main results.