Stochastic models of population growth

1.   Introduction

Natural processes are subject to various random perturbations, so we often use stochastic methods to build their models. In this article we restrict ourselves to single population growth models. We show how a stochastic perturbation is introduced into a deterministic model. We start with the classic situation when we add a diffusion-type noise to a population growth equation. Such models are already quite well studied ^[1,2,3,4,5], but for the sake of completeness we present their important properties in a concise way.

The next class of models is obtained by replacing one or more parameters in deterministic models with stochastic processes ^[6,7,8,9,10]. This case is a bit more difficult, because we generally get degenerate diffusion processes. The last class of models considered here are models with disasters ^{[11,12,13,14]}. The growth of population is described by a deterministic flow with random jumps.

Our aim is to present some mathematical tools and methods for investigating the models under consideration, especially long-time behavior of sample paths and distributions of considered stochastic processes. We study when the distribution of the population size stabilizes or when the population dies out (the trajectories converge to zero or the distributions of the process converge weakly to Dirac's delta at zero) or has unlimited growth.

The considered processes generate semigroups of operators on $L^1$ space, describing the evolution of the densities of their distributions. We are particularly interested in the existence of invariant (stationary) densities and their asymptotic stability. The ergodicity of processes follows from this property, and thus, using the Monte Carlo method, we can easily determine invariant densities using computer simulations.

The organization of the paper is as follows. Section 2 is devoted to population growth models described by one-dimensional stochastic equations of the Itô type. We consider models with environmental and demographic noise. We recall the properties of such models, in particular the long-time behavior of the sample paths and distributions of processes. In this section, we also introduce the concept of partially integral stochastic semigroups and present theorems on their long-time behavior. In particular, we show when the asymptotic stability follows from the existence of an invariant density (Theorem 4) and give conditions under which Foguel's alternative holds (Theorem 5), i.e., the semigroup is either asymptotically stable or is "sweeping" from compact sets. We also use these theorems in later sections of the paper.

In Section 3 we study a model obtained from an ordinary differential equation by replacing the per capita birth rate by a stationary diffusion process. We obtain a degenerate diffusion process $(X_t, Y_t)$ but, using Hörmander's condition, we show that the transition probability function has a density $p$ and, using a method based on support theorems, we prove positivity of $p$ on a sufficiently large set. Then we give conditions when the densities of the process $(X_t, Y_t)$ converge to a stationary density. Such a property is usually proved by checking the tightness of the distributions of the considered process and showing that the stationary measures on the boundary of the domain are repulsive ^[15,16,17]. These conditions can be proved using an appropriately chosen Lyapunov function (Khasminskii function) ^[18]. In our case, the construction of the Lyapunov function is quite difficult, so we use a different technique. First we show that the stochastic semigroup generated by the process $(X_t, Y_t)$ satisfies the Foguel alternative. Then, using an ergodic theorem, we prove the existence of a compact set $C$ satisfying the condition $\text{P}((X_t, Y_t)\in C) \ge \varepsilon > 0$, which guarantees asymptotic stability. We also present numerical simulations of the theoretical results.

In Section 4 we consider a model of population growth with random disasters. This model is obtained by adding sudden decreases in its size at random times to a differential equation describing the growth of the population. The population size distribution satisfies a partial differential equation with a non-local operator describing the jump process. This equation generates a stochastic semigroup, which we investigate using the Foguel alternative. We analyze the long-time behavior of the semigroup in the case where the sudden decrease in population size is proportional to the population size.

2.   Models described by the Itô stochastic differential equations

In this section we recall properties of models given by the Itô stochastic differential equations. Such models are obtained from deterministic models by adding a stochastic perturbation. A deterministic model of a single population is of the form:

where $x(t)$ is the size of the population at time $t$, and $b(x)$ and $l(x)$ are the per capita birth and death rates, respectively. We assume that $b$ and $l$ are $C^1$ nonnegative functions and $b-l$ is a bounded above function. The function $g(x) = (b(x)-l(x))x$ is called the growth rate. Thus we can write Eq (2.1) in the form $\dot x(t) = g(x(t))$. By adding noise to the last equation we obtain the following Itô stochastic equation:

Here and throughout the paper, $W_t$ denotes the standard Wiener process and $\sigma(x)$ is the diffusion coefficient. We usually consider two forms of the function $\sigma(x)$. If we consider environmental noise, then $\sigma(x) = \sigma x$. In the case of demographic noise we assume that $\sigma(x) = \sigma \sqrt{x}$.

In order to study properties of sample paths of the process $(X_t)$ we need the following result being a consequence of the Feynman–Kac formula.

Theorem 1. We assume that $g$ and $\sigma$ are $C^1$-functions in the interval $[\alpha, \beta]$ and $\sigma(x) > 0$ for $x\in [\alpha, \beta]$. Let $X_0 = x_0\in (\alpha, \beta)$ and $\tau$ be the first exit time from the interval $[\alpha, \beta]$. Then $\text{P}(\tau < \infty) = 1$ and

where $\gamma$ is any point from the interval $[\alpha, \beta]$.

Since the properties of the solutions of Eq (2.2) differ significantly depending on the choice of the function $\sigma$ we consider the both cases separately.

2.1.   Environmental noise

Now assume that $\sigma(x) = \sigma x$, $\sigma > 0$. The long-time behavior of sample paths of the process $(X_t)$ depends on whether the following expressions:

are finite. It is worth noting that the expression $J_1(x)$ can be infinite in the case when the function under the integral sign is unbounded at zero (see Example 2). The following results can be obtained from Theorem 1 by the limit passages with $\alpha\to 0$ and $\beta\to\infty$. Let $X_0 = x$, $x\in (0, \infty)$. Then

a) if $J_1(x) = +\infty$ and $J_2(x) = +\infty$, then

b) if $J_1(x) < +\infty$ and $J_2(x) < +\infty$, then

c) if $J_1(x) < +\infty$ and $J_2(x) = +\infty$, then

d) if $J_1(x) = +\infty$ and $J_2(x) < +\infty$, then

Example 2. If $g(x) = g x$, then

Consequently, $J_1(x) < \infty$ if and only if $2g < \sigma^2$ and $J_2(x) < \infty$ if and only if $2g > \sigma^2$. Thus, if $2g = \sigma^2$, then a) holds; if $2g < \sigma^2$, then c) holds; and if $2g > \sigma^2$, then d) holds. These properties also follow from the formula

Moreover, from (2.8) it follows that $\lim_{t\to\infty}\text{P} (X_t\le M) = 1$ if $2g < \sigma^2$; $\lim_{t\to\infty}\text{P} (X_t\ge M) = 1$ if $2g > \sigma^2$; $\lim_{t\to\infty}\text{P} (X_t\le M) = 1/2$ and $\lim_{t\to\infty}\text{P} (X_t\ge M) = 1/2$ if $2g = \sigma^2$ for each $M > 0$.

Theorem 3. We assume additionally that $g$ is a $C^2$-function and $g(x)/x\le \sigma^2/2$ for sufficiently large $x$. Let $\bar g = g'(0)-\sigma^2/2$.

(i) If $\, \bar g < 0$, then $\lim_{t\to\infty} X_t = 0$ a.e.

(ii) If $\, \bar g\ge 0$, then $\liminf_{t\to\infty} X_t = 0$ and $\limsup_{t\to\infty} X_t = \infty$ a.e.

The abbreviation a.e. stands for "almost everywhere".

Proof. From the inequality $g(x)/x\le \sigma^2/2$ it follows that $J_2(x) = +\infty$. We check when $J_1(x)$ is finite. Since $g(0) = 0$, from Taylor's formula it follows that for each $r > 0$ there is $\theta_r\in (0, r)$ such that $2g(r) = 2g'(0)r+g''(\theta_r)r^2$. Thus

where $z_s$ has values in the interval $(0, \max(1, x))$. Since

$J_1(x)$ is finite if and only if $\bar g < 0$. Thus if $\bar g < 0$, then the case c) holds, and if $\bar g\ge 0$, then the case a) holds. □

We now examine the problem of the asymptotic behavior of the distribution of the process $(X_t)$. The variables $X_t$ have densities $u(t, x)$ for $t > 0$, which satisfy the Fokker–Planck equation:

If $X_0$ has the density $u_0$ and $u(t, x)$ is the solution of Eq (2.9), then the formula $P(t)u_0 = u(t, x)$ defines a stochastic semigroup on the space $L^1[0, \infty)$.

We recall that if $(X, \Sigma, m)$ is a $\sigma$-finite measure space, then a family $\{P(t)\}_{t\ge 0}$ of linear operators on $L^1 = L^1(X, \Sigma, m)$ is called a stochastic semigroup if each operator $P(t)$ is positive (if $f\ge 0$ then $P(t)f\ge 0$) and preserves the integral ($\int_X P(t)f(x)\, m(dx) = \int_X f(x)\, m(dx)$ for $f\in L^1$), and the function $t\mapsto P(t)f$ is continuous for each $f\in L^1$.

The semigroup $\{P(t)\}_{t\ge 0}$ is called asymptotically stable if there exists a density $f_*$ such that

for each density $f$. Here $\|\cdot\|$ is the norm in $L^1$, i.e., $\|f\| = \int_X |f(x)|\, m(dx)$. We recall that a measurable function $f\colon X\to \mathbb R$ is a density if $f\ge 0$ and $\|f\| = 1$. From (2.10) it follows immediately that $f_*$ is an invariant density, i.e., $P(t)f_* = f_*$ for each $t\ge 0$.

In order to formulate a theorem on asymptotic stability of stochastic semigroups we need to introduce an auxiliary notion.

A stochastic semigroup $\{P(t)\}_{t\ge 0}$ is called partially integral if there exists a measurable function $k\colon (0, \infty)\times X\times X\to[0, \infty]$, called a kernel, such that

for every density $f$ and

for some $t > 0$.

Theorem 4 (^[19]). Let $\{P(t)\}_{t\ge 0}$ be a partially integral stochastic semigroup. Assume that the semigroup $\{P(t)\}_{t\ge 0}$ has a unique invariant density $f_*$. If $f_* > 0$ a.e., then the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable.

The second important notion which describes the long-time behavior of stochastic semigroups is sweeping. A stochastic semigroup $\{P(t)\}_{t\ge 0}$ is called sweeping with respect to a set $A\in\Sigma$ if for every $f\in L^1$ we have

Further we assume that $(X, \rho)$ is a separable metric space and $\Sigma = \mathcal B(X)$ is the $\sigma$-algebra of Borel subsets of $X$. We consider a partially integral semigroup $\{P(t)\}_{t\ge0}$ which satisfies the following condition:

(K) For every $x_0\in X$ there exists an $\varepsilon > 0$, a $t > 0$, and a measurable function $\eta\ge 0$ such that $\int \eta(x)\, m(dx) > 0$ and

where $B(x_0, \varepsilon) = \{x\in X:\, \, \rho(x, x_0) < \varepsilon\}$ and $\mathbf 1_A$ denotes the indicator function of the set $A$, i.e., $\mathbf 1_A(x) = 1$ if $x\in A$ and $\mathbf 1_A(x) = 0$ if $x\notin A$.

Theorem 5 (^[20]). If a stochastic semigroup $\{P(t)\}_{t\ge0}$ satisfies $\rm (K)$ and $\int_0^{\infty} P(t)f(x) \, dt > 0$ $\rm a.e.$ for every density $f$, then the Foguel alternative holds, i.e., this semigroup is asymptotically stable or sweeping from compact sets.

The semigroup $\{P(t)\}_{t\ge 0}$ generated by Eq (2.9) is integral and has a positive kernel, which implies that it has at most one invariant density. If an invariant density $f_*$ exists then it satisfies the equation

Hence, we easily determine the invariant density:

as long as $C$ can be chosen such that $f_*$ is a density, which reduces to proving the integrability of $f_*$ for $C > 0$. The integrability of the function $f_*$ depends only on its behavior near $0$ and $+\infty$. The study of this integrability can be carried out similarly to the study of the sample paths properties. According to Theorem 4, we have the following result.

Theorem 6. Let $g$ be a $C^2$-function, $\limsup_{x\to\infty} g(x)/x < \sigma^2/2$, and $\, \bar g > 0$. Then the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable.

Example 7 (Logistic model). Now we consider a stochastic version of the logistic model. Let $g(x) = rx(K-x)$, $r > 0$, and $K > 0$. Then $\limsup_{x\to\infty} g(x)/x = -\infty$ and $\bar g = rK-\sigma^2/2$. Thus, if $2rK > \sigma^2$, then the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable. Using the formula (2.14) we obtain that $f_*$ is the gamma density

where $\lambda = \frac{2rK}{\sigma^2}-1$, and $\theta = \frac{\sigma^2}{2r}$. Here $\text{E }X = \lambda\theta = K-\frac{\sigma^2}{2r}$.

Example 8 (Gompertz model). In the study of cell populations (e.g., tumours), we often use the Gompertz model to describe population growth. In such a model, the size of the population satisfies the equation $\dot x = g(x)$, in which $g(x) = rx\ln(K/x)$, $r, K > 0$. Thus a stochastic version of the Gompertz model is the following:

Note that the function $g$ is not differentiable at $x = 0$, so we cannot apply the theory presented earlier. After substituting $Z_t = \ln X_t$ and using the Itô formula we obtain

where

Then

If $V_t = Z_t-c$ and $c = \ln K-\sigma^2/(2r)$, then $V_t$ satisfies the Langevin equation

One of the solutions of Eq (2.17) is the stationary Ornstein–Uhlenbeck process which has the invariant density $v_*(x) = \sqrt{r/\pi\sigma^2}\exp(-rx^2/\sigma^2)$. Since $X_t = \exp(V_t+c)$, the process $(X_t)$ has the log-normal invariant density

Thus, the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable.

Since for any initial distribution of $X_0$ and any $t > 0$ the random variables $X_t$ have densities $u(t, \cdot)$, we finally obtain

for any distribution of $X_0$. Thus we can use the ergodic theorem, which states that if $\varphi\colon [0, \infty)\to \mathbb R$ is a measurable function and $\varphi f_*$ is an integrable function, then

for an a.e. sample path. In particular, for any measurable set $A \subset [0, \infty)$ the mean sojourn time of almost all trajectories in the set $A$ is $\int_Af_*(x)\, dx$.

From Theorem 3 (i) it follows that if $\, \bar g < 0$ and $g(x)/x\le \sigma^2/2$ for sufficiently large $x$, then $\lim_{t\to\infty} X_t = 0$ a.e., which means that the population dies out. Also from Theorem 3 (i) we obtain

One can check that if $\, \bar g = 0$ and $g(x)/x\le \sigma^2/2-\varepsilon$ for some $\varepsilon > 0$ and sufficiently large $x$, then the condition (2.19) also holds, i.e., the population dies out in the sense of distribution.

Another more realistic definition of extinction can be considered by assuming that a population goes extinct if the number of individuals falls below a certain critical level $x_{\text{crit}} > 0$ and then extinction will occur over a finite period of time. In all cases considered so far, this condition is fulfilled. It follows from the fact that if $g(x)/x\le \sigma^2/2$ for sufficiently large $x$, then almost all sample paths of the process $(X_t)$ visit any neighborhood of zero with a probability of one.

For a population to survive with a positive probability, we need to assume that the function $g$ is growing at a greater rate than $\sigma^2 x/2$.

Example 9. Consider the model with $g(x) = g x$ and $\sigma(x) = \sigma x$. We assume that $X_0 = x_0 > x_\text{crit}$ and $\bar g = g-\sigma^2/2 > 0$. Then

This result follows from Theorem 1. We omit the detailed proof.

2.2.   Demographic noise

In the models we have considered so far, the functions $g$ and $\sigma$ were sufficiently smooth and satisfied the condition $g(0) = \sigma(0) = 0$. As a result, the diffusion process never hits the point $0$. In such a case, we say that the point $0$ belongs to the unattainable boundary for the diffusion process (if the process hits the boundary with positive probability, then this boundary is called attainable). We check that if the noise is demographic, i.e., $\sigma(x) = \sigma \sqrt{x}$, then $0$ is an attainable point.

We start with a particular model:

If $g = 0$, then the process $(X_t)$ is called the Feller diffusion.

We show that the process $(X_t)$ satisfying (2.20) hits $0$ with positive probability. In order to do it consider Eq (2.20) in the interval $[\alpha, \beta]$ with $0 < \alpha < x_0 < \beta < \infty$. According to formula (2.3), the probability that the process $(X_t)$ exits through the right end of the interval $[\alpha, \beta]$ is $\Phi(x_0)/\Phi(\beta)$, where

for $g\ne 0$, and $\Phi(x) = x-\alpha$ for $g = 0$. Hence

On passing to the limit $\beta \to+\infty$ we get the following:

On passing this time to the limit $\alpha\to 0^+$ we obtain

so the probability that the process $(X_t)$ hits zero equals $e^{-2gx_0/\sigma^2}$ when $g > 0$, while it is $1$ when $g\le 0$.

Remark 10. Eq (2.20) has the following interesting property. If $(X_t^1)$ and $(X_t^2)$ are independent processes defined by this equation satisfying the initial conditions $X_0^i = x_i$, $i = 1, 2$, then the process $(X_t)$ satisfying this equation with initial condition $X_0 = x_1+x_2$ has the same distribution as the process $(X_t^1+X_t^2)$. This fact is related to the demographic noise property that the stochastic disturbance of the entire population is the sum of the independent disturbances of the subpopulations.

We now consider the general model:

where $g$ is a $C^1$-function in the interval $[0, \infty)$. It is easy to check that if $\limsup_{x\to\infty} g(x) < \sigma^2/2$, then $\Phi(\infty) = \infty$, which implies that the process $(X_t)$ hits zero with a probability of one. If $\liminf_{x\to\infty} g(x) > \sigma^2/2$ and $g'(0) > 0$, then $\Phi(\infty) < \infty$ and the process $(X_t)$ hits zero with a probability less than one. Precisely, sample paths with positive probability hit zero and with positive probability converge to $+\infty$.

Comparing the properties of models with environmental noise and demographic noise we deduce that for small populations the demographic noise is more dangerous than the environmental noise. The opposite is true for large populations, i.e., the demographic noise has less impact on the population.

3.   A model with a stochastic parameter

A new class of stochastic models can be obtained if we replace a parameter in a deterministic model by a stochastic process. We begin with a simple deterministic model with a constant per capita birth rate $b$:

We assume that the external noise changes the birth rate $b$, and so we replace $b$ by a stationary Markov process $(Y_t)$ which satisfies the Itô equation

Since the death rate is a positive number, we assume that the sample paths of this process are positive functions. In order to keep positive values of sample paths, we assume that $a$ is a $C^1$-function defined on the interval $[0, \infty)$. According to Theorem 6, so that there is a stationary process $(Y_t)$ satisfying Eq (3.2), it is sufficient to assume that $a$ is a $C^2$-function, $a(0) > \sigma^2/2$, and $\limsup_{x\to\infty} a(x) < \sigma^2/2$.

We replace $b$ in Eq (3.1) by the process $(Y_t)$. Then the solution $x(t)$ of this equation is also a stochastic process, denoted by $(X_t)$, and the growth of the population is described now by the following equation:

We also assume that the process $(Y_t)$ has a finite mean value $\text{E }Y$. It is easy to check that if $\limsup_{x\to\infty} a(x) < 0$, then $\text{E }Y < \infty$. We denote by $(\Omega, \mathcal F, \text{P})$ the probability space on which the process $(X_t, Y_t)$ is defined.

Remark 11. One can consider also models with values of sample paths of $(Y_t)$ in some interval $(\alpha, \beta)$, but then we need to consider Eq (3.2) with drift and diffusion coefficients having zeroes at $\alpha$ and $\beta$, and fulfilling an appropriate condition guaranteeing the existence of a stationary process.

We precede our study of the long-time behavior of the solutions of Eq (3.1) with a remark on "medium-time behavior". Let a function $x\colon [0, T]\to (0, \infty)$ be a sample path of the process $(X_t)$, i.e., $x(t) = X_t(\omega)$ for some $\omega\in\Omega$. Since the sample paths of the processes $(X_t)$ and $(Y_t)$ are positive functions, according to Eq (3.3), the function $x(t)$ satisfies the following condition:

Of course not all functions satisfying (3.4) are sample paths of $(X_t)$, but we have the following property. For each $\varepsilon > 0$, $T\ge 0$, and $x$ satisfying (3.4) we have

This property follows immediately from the fact that if $y\colon [0, T]\to (0, \infty)$ is any continuous function and $\delta > 0$, then

Thus, we see that any function satisfying condition (3.4) approximately realizes one of the possible population growth scenarios.

The solution of problem (3.3) depends on the process $(Y_t)$, but can only be expressed by explicit formulae in a few cases. The simplest case is when $l$ is a constant, say $l(x) = \rho$. Then

According to the ergodic theorem $\lim_{t\to\infty}\frac 1t\int_0^tY_s\, ds = \text{E }Y$ a.e. Thus, the a.e. sample path of the process $(X_t)$ approaches exponentially $+\infty$ when $\text{E }Y > \rho$, and $0$ when $\text{E }Y < \rho$. In this model the long-time behavior is practically the same as in the deterministic model with the per capita growth rate $\text{E }Y-\rho$.

Now we consider a logistic model with $l(x) = \rho+\gamma x$. Then

Unfortunately, the application of the ergodic theorem to the process $(Y_t)$ is not sufficient to examine the properties of the sample paths of the process $(X_t)$. Further investigation requires the use of methods related to degenerate diffusion. The following theorem allows us to study the properties of distributions and sample paths of processes $(X_t, Y_t)$ for a rather broad class of models including the logistic one. We use the notation $\mathbb R_+^2 = [0, \infty)^2$ and $\mathbb R_+^{2, \circ} = (0, \infty)^2$.

Theorem 12. Let $a$ and $l$ be $C^{\infty}$-functions such that

Then for each $t > 0$ the distribution of $(X_t, Y_t)$ has a density $f_t$ and there exists a density $f_*$ such that $\lim_{t\to\infty}\|f_t-f_*\| = 0$ in $L^1(\mathbb R_+^{2, \circ})$.

We explain the assumptions of Theorem 12. Inequalities $a(0) > \sigma^2/2$ and $\limsup_{x\to\infty} a(x) < 0$ imply that the process $(Y_t)$ is stationary and has a finite mean value. Inequalities containing $\text{E }Y$ imply that there exists a compact subset $C$ of $\mathbb R_+^{2, \circ}$ such that $\text{P}((X_t, Y_t)\in C) > \varepsilon > 0$.

The plan of the proof is the following. First, we prove that the transition probability function $\mathcal P(t, (x_0, y_0), A)$ of the process $(X_t, Y_t)$ has a continuous density $p(t, (x_0, y_0), (x, y))$ for each $t > 0$. Next we check that for arbitrary positive numbers $x_0, y_0, x_1, y_1$ there exists $T > 0$ such that $p(T, (x_0, y_0), (x_1, y_1)) > 0$. Then the stochastic semigroup $\{P(t)\}_{t\ge 0}$ on the space $L^1(\mathbb R_+^{2, \circ})$ generated by the process $(X_t, Y_t)$ satisfies the assumptions of Theorem 5. Thus the Foguel alternative holds. Since we find a compact set which does not satisfy sweeping condition (2.11), the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable.

Lemma 13. The transition probability function $\mathcal P(t, (x_0, y_0), A)$ has a density $p(t, (x_0, y_0), (x, y))$ and $p\in C^\infty((0, \infty)\times \mathbb R_+^{2, \circ} \times \mathbb R_+^{2, \circ})$.

Proof. First we convert the Itô system (3.2)–(3.3) into the Stratonovitch one ^[21,22]:

We prove the lemma by checking Hörmander's condition (cf. ^[23]). Let $\sigma_0(x, y) = ((y-l(x))x, (a(y)-\frac12\sigma^2)y)$ and $\sigma_1(x, y) = (0, \sigma y)$. It is enough to check that the vector $\sigma_1(x, y)$ and the Lie bracket $[\sigma_0, \sigma_1](x, y)$ span the space $\mathbb R^2$ for each $(x, y)\in \mathbb R_+^{2, \circ}$. We recall that

We have

and hence $[\sigma_0, \sigma_1](x, y) = -\sigma (xy, a'(y)y^2)$. It is clear that vectors $\sigma_1(x, y)$ and $[\sigma_0, \sigma_1](x, y)$ are linearly independent for each $(x, y)\in \mathbb R_+^{2, \circ}$. □

There is one important difference between a non-degenerate diffusion and a degenerate diffusion which satisfies Hörmander's condition, namely, the kernel $p$ is strictly positive if diffusion is non-degenerate, but in the degenerate case the kernel $p$ can vanish on some subsets. We check where the kernel $p$ is positive using a method based on support theorems (^[24,25,26]):

Theorem 14. Fix a point $(x_0, y_0)\in \mathbb R_+^{2, \circ}$ and assume that we can find a function $\phi\in L^2([0, T];\mathbb R)$ such that there exists a solution of the system

If $(x, y) = (x_{\phi}(T), y_{\phi}(T))$ and the Frechét derivative of the function $\Psi(h) = (x_{\phi+h}(T), y_{\phi+h}(T))$ in $h_0\equiv 0$ has the rank $2$, then $p(T, (x_0, y_0), (x, y)) > 0$.

It is worth noting that this is a strong version of the support theorem for a degenerate diffusion. The topological support of the measure $\mathcal P(T, (x_0, y_0), \cdot)$ coincides with the closure of the set of all points $(x, y)$ which can be connected to $(x_0, y_0)$ at time $T$ by means of solving system (3.9) with an appropriately chosen function $\phi$.

Lemma 15. For any quadruple $(x_0, y_0, x_1, y_1)$ of positive numbers there exists $T > 0$ such that $p(T, (x_0, y_0), (x_1, y_1)) > 0$.

Proof. Fix positive numbers $x_0, x_1, y_0, y_1$. We show that there exist $T > 0$ and a continuous function $\phi\colon [0, T]\to\mathbb R$ such that the system

has a solution satisfying the following boundary conditions:

Observe that it is sufficient to find $T > 0$, a positive function $x\in C^1[0, T]$ satisfying the inequality

the boundary conditions $x(0) = x_0$, $x(T) = x_1$ and the consistency conditions $\dot x(0) = (y_0-l(x_0))x_0$, $\dot x(T) = (y_1-l(x_1))x_1$. Given a function $y$, we find the function $\phi$ by solving Eq (3.11). Let $F\colon (0, \infty)\to \mathbb R$ be a function given by

Let $z(t) = F(x(t))$. We can write inequality (3.13) as $\dot z(t) > -1$ for $t\in [0, T]$. We have $x(0) = x_0 \Leftrightarrow z(0) = 0$, $x(T) = x_1 \Leftrightarrow z(T) = F(x_1)$ and the consistency conditions can be written as

Since $\dot z(0) > -1$ and $\dot z(T) > -1$, it is not difficult to find $T > 0$ and a $C^1$-function $z\colon [0, T]\to\mathbb R$ which fulfils all the required conditions. Since $F$ is a smooth strictly increasing function we can define $x\colon [0, T]\to\mathbb R$ as $x(t) = F^{-1}(z(t))$ for $t\in [0, T]$.

The Frechét derivative of the function $\Psi$ can be found by means of the perturbation method for ordinary differential equations. Let $\mathbf x = (x, y)$, $\mathbf x_{\phi} = (x_{\phi}, y_{\phi})$,

and let $Q(t, t_0)$, for $T\ge t\ge t_0\ge 0$, be a matrix function such that $Q(t_0, t_0) = I$ and $\dfrac{\partial Q(t, t_0)}{\partial t} = \Lambda(t)Q(t, t_0)$. Then

Let $\varepsilon\in(0, T)$ and $h_{\varepsilon} = \mathbf 1_{[T-\varepsilon, T]}$. Since $Q(T, s) = I+\Lambda(T)(T-s)+o(T-s)$, from (3.15) we obtain

We have $\mathbf x_{\phi}(T) = \mathbf x$. Let $\mathbf v_1 = \sigma_1(\mathbf x_{\phi}(T))$ and $\mathbf v_2 = \Lambda(T) \sigma_1(\mathbf x_{\phi}(T))$. Then

Since the vectors $\mathbf v_1$ and $\mathbf v_2$ are linearly independent, the vectors $\Psi'(h_{\varepsilon})$ and $\Psi'(h_{2\varepsilon})$ are also linearly independent. Thus the derivative $\Psi'$ has rank 2 and finally $p(T, (x_0, y_0), (x, y)) > 0$. □

Let $\bar Y_t = Y_t-\text{E }Y$ and $Z_t = \ln X_t$. Then

where $c(z) = l(e^z)-\text{E }Y$. From our assumptions it follows that there exist $\delta > 0$ and $\bar z > 0$ such that $c(z) > \delta$ for $z\ge \bar z$ and $c(z) < -\delta$ for $z\le -\bar z$.

Lemma 16. Assume that the random variable $Z_0$ has values in the interval $[-\bar z, \bar z]$. Then there exist $T_0 > 0$, and $M > \bar z$ such that $\text{P}(\omega\in\Omega \colon |Z_t(\omega)|\le M)\ge 1/2$ for $t\ge T_0$.

Proof. The process $(\bar Y_{-t})$ is also stationary and $\text{E }\bar Y = 0$. From the ergodic theorem applying to the process $(\bar Y_{-s})$, $s\ge 0$, it follows that

Since a.e. convergence implies convergence in probability, there exists $T_0 > 0$ such that

Now let $t\ge T_0$. Then from the last formula it follows that

Since the process $(\bar Y_t)$ is stationary and its sample paths are continuous functions, there exists a constant $K > \delta$ independent of $t$ such that

Let $\Omega_0$ be a subset of $\Omega$ containing all elementary events $\omega$ satisfying both conditions (3.18) and (3.19). Then $\text{P}(\Omega_0)\ge 1/2$. Define $M = \bar z+(K-\delta)T_0$. We check that $|Z_t(\omega)| < M$ for $\omega\in \Omega_0$. Let $\omega\in \Omega_0$ and $\tau(\omega)\le t$ be such that $Z_{\tau(\omega)}\in [-\bar z, \bar z]$ and $Z_{s}\notin [-\bar z, \bar z]$ for $s\in (\tau(\omega), t]$. Without loss of generality we can assume that $\tau(\omega) < t$. If $Z_{\tau(\omega)} = \bar z$, then

Let $I_1 = \int_{\tau(\omega)}^t \bar Y_s(\omega)\, ds$. Then from (3.18) and (3.19) we obtain

Since $c(Z_s(\omega))\ge \delta$ for $s\in [\tau(\omega), t]$, we have $\int_{\tau(\omega)}^t c(Z_s(\omega))\, ds\ge \delta(t-\tau(\omega))$. Thus $Z_t(\omega)\le \bar z+(K-\delta)T_0 = M$ for $\omega\in\Omega_0$. If $Z_{\tau(\omega)} = -\bar z$, then we analogously obtain that $Z_t(\omega)\ge -M$ for $\omega\in\Omega_0$. □

Proof of Theorem 12. First we check that the stochastic semigroup $\{P(t)\}_{t\ge 0}$ on the space $L^1(\mathbb R_+^{2, \circ})$ generated by the process $(X_t, Y_t)$ satisfies the assumptions of Theorem 5. The semigroup $\{P(t)\}_{t\ge 0}$ is given by the formula

According to Lemma 15 for any quadruple $(x_0, y_0, x, y)$ of positive numbers there exists $t > 0$ such that $p(t, (x_0, y_0), (x, y)) > 0$. This and the continuity of $p$ imply condition (K). Also Lemma 15 and formula (3.20) imply that $\int_0^{\infty}P(t)f(x, y)\, dt > 0$ for every density $f$ and for all $(x, y)\in \mathbb R_+^{2, \circ}$. Thus the Foguel alternative holds.

Let $f_1(x)$ be a density supported on the interval $[e^{-\bar z}, e^{\bar z}]$. If $f_1$ is the density of the distribution of $X_0$, then according to Lemma 16 we have $\text{P}\big(X_t\in \big[e^{-M}, e^{M}\big]\big)\ge 1/2$ for $t\ge T_0$. Let $f_*(y)$ be the invariant density for the process $(Y_t)$. Then there exist positive numbers $y_1$ and $y_2 > y_1$ such that $\text{P}\big(y_1\le Y_t\le y_2\big)\ge 3/4$ for each $t\ge 0$. Let $f(x, y) = f_1(x)f_*(y)$ and $C = \big[e^{-M}, e^{M}\big]\times \big[y_1, y_2\big]$. Then

Thus the semigroup $\{P(t)\}_{t\ge 0}$ is not sweeping from $C$ and, consequently, the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable. □

We now present the density plots of the joint stationary density of the process $(X_t, Y_t)$ and the marginal densities of the distributions of the processes $(X_t)$ and $(Y_t)$ when $l(x) = c+x$, $a(y) = 5-y$, and $\sigma = 1$. Since the process $(Y_t)$ is a stationary solution of the logistic model with $g(y) = y(5-y)$ it has the invariant density $f(y) = \frac{2^9}{\Gamma(9)} y^{8}e^{-2y}$ and $\text{E }Y = 4.5$ (see Figure 1).

The plots of the joint stationary distribution of the process $(X_t, Y_t)$ and the marginal density of $(X_t)$ are given in Figure 2 for $c = 1$ and in Figure 3 for $c = 4$. Observe that since $\ln X_t-\ln X_0 = Y_t-(c+X_t)$ we have $\text{E }Y = c+\text{E }X$ and consequently $\text{E }X = 3.5$ in the case of $c = 1$ and $\ E X_t = 0.5$ in the case of $c = 4$.

The graphs in Figures 2 and 3 were prepared using MATLAB. Each plot was obtained from $1000$ simulations of sample paths up to time $T = 100$ with the time step size $\Delta t = 0.001$.

4.   Population growth with random disasters

We now assume that a population grows according to the equation

but from time to time, the population encounters a catastrophe that causes its size to drop from $x$ to $y$, where $y$ is distributed according to a transition probability function $\mathcal P(x, A)$ (see Figure 4). We assume that the measure $\mathcal P(x, \cdot)$ has a density $p(x, y)$, i.e., $p\ge 0$, $\int_0^x p(x, y)\, dy = 1$, and $\mathcal P(x, A) = \int_A p(x, y)\, dy$. It is also assumed that the moments of disasters are distributed according to a Poisson process $(N_t)_{t\ge 0}$ with rate $\Lambda > 0$.

We assume that $g$ is a $C^1$-function, $g(0) = 0$, and $g(x) > 0$ for $x\in (0, \alpha)$, where $\alpha\in (0, \infty]$. If $\alpha < \infty$, then we assume that $g(\alpha) = 0$. We also assume that $\int_{\varepsilon}^{\alpha}\frac{dx}{g(x)} = \infty$ for $\varepsilon > 0$, which guarantees the existence of the solutions of Eq (4.1) for all $t > 0$. We assume that $p(x, y)$ is a continuous function for $0 < y < x < \alpha$.

The growth of the population is described by a piecewise deterministic Markov process $(X_t)$, $t\ge 0$, with random jumps at the moments of disasters. Assume that the distribution of random variable $X_0$ has a density $f$. Then for all $t > 0$ the distribution of $X_t$ has also a density $u(t, \cdot)$ and $u$ satisfies the following equation:

where $Q$ is a stochastic operator on the space $L^1: = L^1(0, \alpha)$ given by $Qf(y) = \int_0^{\alpha} p(x, y)f(x)\, dx$.

Equation (4.2) can be written as an evolution equation $\dot u(t) = \mathcal Au(t)$ in the space $L^1$, where $\mathcal A = \mathcal A_0+\Lambda Q-\Lambda I$, the operator $\mathcal A_0$ is given by

and the domain of both operators $\mathcal A$ and $\mathcal A_0$ is the set $\{f\in L^1\colon \mathcal A_0f\in L^1\}$. We denote by $\pi_tx_0$ the solution $x(t)$ of Eq (4.1) with the initial condition $x(0) = x_0\ge 0$. The operator $\mathcal A_0$ generates a stochastic semigroup $\{P_0(t)\}_{t\ge 0}$ on the space $L^1$ given by the formula

From the Dyson–Phillips expansion the operator $\mathcal A$ also generates a stochastic semigroup $\{P(t)\}_{t\ge 0}$ on the space $L^1$ given by the formula

where

It is easy to check that $P_1(t)$ is an integral stochastic operator with a continuous kernel for $t > 0$, which implies that the semigroup $\{P(t)\}_{t\ge 0}$ satisfies condition (K). Moreover, the semigroup $\{P(t)\}_{t\ge 0}$ satisfies condition $\int_0^{\infty}P(t)f(x)\, dt > 0$ for all $x\in (0, \alpha)$ and each density $f$. The proof of this result is a little bit technical and we omit it here, but it is a consequence of the fact that the process $(X_t)$ joins an arbitrary two points from the interval $(0, \alpha)$. According to Theorem 5 we have:

Corollary 17. The semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable or sweeping from compact subsets of the interval $(0, \alpha)$.

Under the assumptions made so far, the following cases are possible:

(a) the semigroup is asymptotically stable,

(b) the semigroup is sweeping to zero, i.e., the distributions of the process $(X_t)$ are weakly converging to $\delta_0$ (population extinction),

(c) the semigroup is sweeping to $\alpha$, i.e., the distributions of the process $(X_t)$ are weakly convergent to $\delta_{\alpha}$ when $\alpha < \infty$ or $\lim_{t\to\infty} \text{P}(X_t\ge M) = 1$ for all $M > 0$ when $\alpha = \infty$,

(d) the distribution of the process $(X_t)$ is in part swept to zero and in part swept to $\alpha$.

One method of investigating which of these cases occurs is to examine the stationary solutions of Eq (4.2). Observe that a stationary solution $f$ satisfies the equation

This equation can be written in the form

Then $\Lambda(\Lambda I-\mathcal A_0)^{-1}$ is also a stochastic operator given by the formula

Thus $f$ is a solution of the following integral equation:

where

If a density $f_*$ is a solution of Eq (4.5), then the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable. If Eq (4.5) has a positive but not integrable solution $f_*$, then one can prove that if $f$ is integrable in a neighborhood of zero, then condition (c) holds, and if $f$ is integrable in a neighborhood of $\alpha$, then condition (b) holds.

Another method of checking which condition holds is to use the Khasminskii function (also known as the Lyapunov function) (see e.g., [28, page 128] for a general result). The process $(X_t)$ has the infinitesimal generator

where $Q^*V(x) = \int_0^x p(x, y)V(y)\, dy$. Assume that there exist a $C^1$-function $V\colon [0, \alpha)\to [0, \infty)$, constants $\varepsilon, M > 0$, and a compact set $C\subset (0, \alpha)$ such that

Then the semigroup $\{P(t)\}_{t\ge 0}$ is not sweeping from the set $C$ and, consequently, this semigroup is asymptotically stable. Using this criterion, we prove the following theorem.

Theorem 18. Assume that $\alpha = \infty$ and that if $x$ is the size of the population before the disaster, then the population size after the disaster is $x\eta$, where the random variable $\eta$ has values in the interval $[0, 1]$ and has a density $q\colon [0, 1]\to [0, \infty)$. Assume that the limit $g'(\infty) = \lim_{x\to\infty} g(x)/x$ exists. If

then the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable.

Proof. We have $p(x, y) = \frac 1xq\left(\frac yx\right)$ for $0\le y\le x$ and $p(x, y) = 0$ otherwise. Let $V\colon (0, \infty)\to (0, \infty)$ be a $C^1$-function such that $V(x) = x^{-\beta}$ for $x$ sufficiently close to zero, and $V(x) = x^{\gamma}$ for sufficiently large $x$, where $\beta$ and $\gamma$ are positive numbers. Then

for $x$ sufficiently close to zero and

for sufficiently large $x$. Since

we find $\beta$ and $\gamma$ sufficiently close to zero and $\varepsilon > 0$ such that $\mathcal LV(x)\le -\varepsilon$ outside some compact set $C\subset (0, \infty)$. This implies that the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable. □

Example 19. Assume that $q(x) = (\alpha+1)x^{\alpha}\mathbf 1_{[0, 1]}(x)$ and $\alpha\ge 0$. Then $\int_0^1 q(r)\ln r\, dr = -(1+\alpha)^{-1}$. From (4.8) it follows that if

then the semigroup $\{P(t)\}_{t\ge 0}$ is asymptotically stable, which means that the distribution of the population size tends toward an invariant distribution.

Remark 20. One can check that if

the population goes extinct and if

the size of the population grows to $\infty$.

In the next theorem, we examine more closely the long-time behavior of the population size distribution when $g(x) = cx$ and one of the conditions in Remark 20 is satisfied.

Theorem 21. Let $g(x) = cx$, $c > 0$, and $p(x, y) = \frac 1xq\left(\frac yx\right)\mathbf 1_{[0, x]}(y)$ for a density $q\colon [0, 1]\to[0, \infty)$. Assume that $m = \int_0^1q(x)\ln x\, dx < \infty$ and $\sigma^2 = \int_0^1q(x)\ln^2 x\, dx < \infty$. If $F_t$ is the cumulative distribution function (CDF) of $X_t$, then

where $\Phi$ is the CDF of the standard normal distribution $\mathcal N(0, 1)$ and the convergence in (4.9) is uniform in $x$.

Proof. First observe that the operator $Q$ and the semigroup $\{P_0(t)\}_{t\ge 0}$ commute, i.e., $QP_0(t) = P_0(t)Q$ for all $t\ge 0$. Indeed, since $P_0(t)f(x) = e^{-ct}f(e^{-ct}x)$ and $Qf(x) = \int_x^{\infty} y^{-1}q(x/y)f(y)\, dy$, we have

We check by induction that the operators $P_n(t)$ and $Q$ also commute and that

for all $n\ge 0$ and $t\ge 0$. From (4.3) we obtain

Let $\xi_0, \eta_1, \eta_2, \dots$ be a sequence of independent random variables such that $\xi_0$ has the density $f$ and each random variable $\eta_i$ has the density $q$. Then the random variable $\xi_n = \xi_0\eta_1\cdots\eta_n$ has the density $Q^nf$. Let $(N_t)_{t\ge 0}$ be the Poisson process independent from random variables $\xi_0, \eta_1, \eta_2, \dots$ and with rate $\Lambda$. Then the random variable $X_t = e^{ct}\xi_{N_t}$ has the density $P(t)f$. Let $Y_n = \ln \eta_n$ for $n\ge 1$ and $Z_t = \ln X_t$. Then

and $(Z_t-ct)_{t\ge 0}$ is a compound Poisson process. We have $m = \text{E }Y_n$ and $\sigma^2 = {\rm Var\, } Y_n+m^2$. Then from the central limit theorem for the compound Poisson process

uniformly in $x$. Since $Z_t = \ln X_t$, from (4.11) we obtain (4.9). □

Thus the long-time behavior of $\int_0^x P(t)f(y)\, dy$ is the same as the CDF of the log-normal distribution ${\rm LN}((c+\Lambda m)t, \sigma^2\Lambda t)$. From (4.9) it follows that for each $M > 0$ we have $\lim_{t\to\infty}\text{P} (X_t\ge M) = 1$ when $c+\Lambda m > 0$, and $\lim_{t\to\infty}\text{P} (X_t\le M) = 1$ when $c+\Lambda m < 0$, i.e., the population dies out in the sense of distribution.

Example 22. If $q(x) = (\alpha+1)x^{\alpha}\mathbf 1_{[0, 1]}(x)$, $\alpha\ge 0$, then $m = -(1+\alpha)^{-1}$. Consequently, if $\Lambda > c(1+\alpha)$, then the population dies out in the sense of distribution and if $\Lambda < c(1+\alpha)$ the population has unlimited growth.

5.   Discussion and conclusions

In the present paper we have investigated three types of stochastic models of a single population growth. Models with diffusion-type environmental and demographic noise are studied in Section 2. We study long-time behavior of sample paths and densities of distribution and give conditions under which the population will survive and when it will become extinct. Comparing the properties of models with environmental noise and demographic noise we deduce that for small populations the demographic noise is more dangerous than the environmental noise. The opposite is true for large populations, i.e., the demographic noise has less impact on the population. In the case of demographic noise, the population dies out with positive probability for each growth function.

We also investigate when the process describing population size has an asymptotically stable invariant density using simple and convenient tools related to partially integral semigroups of stochastic operators (see Theorems 2.4 and 2.5). These methods can also be applied to more advanced models, for example those considered in Sections 3 and 4.

In Section 3 we studied the model obtained from the deterministic one by replacing the birth rate with a mean-reverting stochastic process ^[6]. We gave sufficient conditions for the existence and asymptotic stability of invariant density. The proof of the main result of Theorem 12 is quite advanced and required the development of new tools that differ from the usual methods for proving asymptotic stability based on the construction of Lyapunov functions. We also ran computer simulations that show how the density of the process changes depending on the value of the death rate at zero.

A model of population growth with random disasters was considered in Section 4. The results obtained make it possible to decide when the population size stabilizes, dies out, or has unlimited growth. It should be noted that we can consider more general models of population growth with random disasters. For example Eq (4.1) can be replaced by a stochastic equation (2.2) and the transition probability function $\mathcal P(x, A)$ can be arbitrary. Then $(X_t)$ is still a Markov process and it can be studied by means of the methods already presented ^[27,28]. More advanced models are those in which the moments of disasters are not distributed according to a Poisson process. For example we can consider a model in which the time between successive disasters is an arbitrary random variable whose distribution depends on the size of the population after the previous disaster. Since $(X_t)$ is now no longer a Markov process, the study of its asymptotic properties requires the use of slightly different methods ^[29].

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare there is no conflict of interest.