Quasistationarity in populations that are subject to large-scale mortality or emigration

by

Phil Pollett

Department of Mathematics The University of Queensland

WHAT ARE WE MODELLING?

Populations which are subject to crashes. Dramatic losses can occur due to

- disease (eg a new virus)

- food shortages (eg overbrowsing)

- significant changes in climate.

Quasistationary behaviour. Such populations can survive for long periods before extinction occurs and can settle down to an apparently stationary regime.

Our goal. We seek to model this behaviour in order to properly manage these populations: to predict persistence times and to estimate population size.

Our model. The birth-death and catastrophe process predicts eventual extinction, but the time till extinction can be very long. The stationarity exhibited by these populations over any reasonable time scale can be explained using a quasistationary distribution.

THE MODEL

We use a continuous-time Markov process (X(t), t ≥ 0), where X(t) is the population size at time t, with transition rates (q_{jk}, j, k ≥ 0) given by

q_{j,j+1} = jρa,  j ≥ 0,
q_{j,j} = −jρ,  j ≥ 0,
q_{j,j−i} = jρb_{i},  j ≥ 2, 1 ≤ i < j,
q_{j,0} = jρ_{∑i≥j} b_{i},  j ≥ 1,

with the other transition rates equal to 0. Here, ρ > 0, a > 0, b_{i} > 0 for at least one i in C = {1, 2, … }, and a + ∑_{i≥1} b_{i} = 1.

Interpretation. For j ≠ k, q_{jk} is the instantaneous rate at which the population size changes from j to k, ρ is the per capita rate of change and, given a change occurs, a is the probability that this results in a birth and b_{i} is the probability that this results in a catastrophe of size i (corresponding to the death or emigration of individuals).

SOME PROPERTIES

The state space. Clearly 0 is an absorbing state (corresponding to population extinction) and C is an irreducible class.

Extinction probabilities. If α_{i} is the probability of extinction starting with i individuals, then α_{i} = 1 for all i ∈ C if and only if D (the expected increment size), given by

D = a −_{∑i≥1} ib_{i} = 1 −_{∑i≥1} (i + 1)b_{i},

is less than 0 (the subcritical case) or equal to 0 (the critical case).

In the supercritical case (D > 0), the extinction probabilities can be expressed in terms of the probability generating function (pgf)

f(s) = a +_{∑i≥1} b_{i}s^{i+1} , |s| < 1.

We find that

∑i≥1 α_{i}s^{i} = s/(1 − s) − Ds/b(s),

where b(s) = f(s) − s.

QUASISTATIONARY

DISTRIBUTIONS

In order to describe the long-term behaviour of the process, we shall use two types of quasistationary distribution (QSD), called Type I and Type II, corresponding to the limits:

lim_{t→∞} Pr(X(t) = j|X(0) = i, X(t) > 0,

X(t + r) = 0 for some r > 0),

lim_{t→∞} lim_{s→∞} Pr(X(t) = j|X(0) = i, X(t + s) > 0,

X(t + s + r) = 0 for some r > 0),

where i, j ∈ C. Thus, we seek the limiting probability that the population size is j, given that extinction has not occurred, or (in the second case) will not occur in the distant future, but that eventually it will occur; we have conditioned on eventual extinction to deal with the supercritical case, where this event has probability less than 1.

THE EXISTENCE OF QSDS

Consider the two eigenvector equations

\begin{matrix}_{∑i∈C} m_{i}q_{ij} = −μm_{j}, j ∈ C, \\ _{∑j∈C} q_{ij}x_{j} = −μx_{i}, i ∈ C, \end{matrix}

where μ ≥ 0 and C is the irreducible class.

In order that both QSDs exist, it is necessary that these equations have strictly positive solutions for some μ > 0, these being the positive left and right eigenvectors of Q_{C} (the transition-rate matrix restricted to C) corresponding to a strictly negative eigenvalue −μ.

Let λ be the maximum value of μ for which positive eigenvectors exist (λ is known to be finite), and denote the corresponding eigenvectors by m = (m_{j}, j ∈ C) and x = (x_{j}, j ∈ C).

THE EXISTENCE OF QSDS

Proposition 5.1 of [1]^{∗} can be restated for our purposes as follows (all sums are over k in C):

Proposition 1. Suppose that Q is regular.

(i) If ∑ m_{k}x_{k} converges, and either ∑ m_{k} converges or {x_{k}} is bounded, then the Type II QSD exists and defines a proper probability distribution π^{(2)} = (π_{j}^{(2)} , j ∈ C) over C, given by

π_{j}^{(2)} = \frac{m_{j}x_{j}}{∑ m_{k}x_{k}}, j ∈ C.

(ii) If in addition ∑ m_{k}α_{k} converges, then the Type I QSD exists and defines a proper probability distribution π^{(1)} = (π_{j}^{(1)} , j ∈ C) over C, given by

π_{j}^{(1)} = \frac{m_{j}α_{j}}{∑ m_{k}α_{k}}, j ∈ C.

∗[1] Pollett, P. (1988) Reversibility, invariance and μinvariance, Adv. Appl. Probab. 20, 600-621.

GEOMETRIC CATASTROPHES

We first examine the important special case b_{j} = b(1 − q)q^{j−1}, j ≥ 1, where b > 0, 0 < q < 1 and a + b = 1. Thus, given a jump occurs, it is a birth with probability a or a catastrophe with probability b, and the size of the catastrophe is determined by a geometric distribution.

We need to solve the right and left eigenvector equations. These are, respectively, for j ≥ 1,

jρax_{j+1} − (jρ − μ)x_{j} + j ; ∑ ; k=0 jρb_{j−k}x_{k} = 0,

and, for k ≥ 1,

(k − 1)ρam_{k−1} − (kρ − μ)m_{k}
+ _{j≥k+1}^{∑} jρb_{j−k}m_{j} = 0,

with the understanding that x_{0} = m_{0} = 0.

We find that D = a − b/(1 − q), and that the maximum value of μ for which there exist strictly positive left and right eigenvectors is

λ = \begin{cases}ρD, & if D > 0, \\ 0, & if D = 0, \\ −ρD(1 − q)a/rb, & if D < 0, \end{cases}

where r = a/(b + qa).

We deduce immediately that no QSD exists in the critical case (D = 0). However, in both the supercritical (D > 0) and subcritical (D < 0) cases, both the Type I and the Type II QSDs exist.

Supercritical case. When D > 0 we find that π_{j}^{(1)} = (1 − 1/r)(1/r)^{j−1} (note that r > 1 since D > 0) and

π_{j}^{(2)} = (1 − q)(r − 1)^{(1+γ)} r^{−(γ+j)}

\binom{γ + j − 1}{j − 1}\frac{γ + (1 − qr)(j − 1)}{γ + j − 1} ,

where γ = (1−qr)/(1−q); note that 0 < γ < 1.

Subcritical case. When D < 0 we find that π_{j}^{(1)} = (1 − r)r^{j−1}, where as above r = a/(b + qa), but now r < 1 since D < 0. We also find that

π_{j}^{(2)} = (1 − qr)(1 − r)^{(1+δ)} r^{j−1}

\binom{δ + j − 1}{j − 1}\frac{δ + (1 − q)(j − 1)}{δ + j − 1} ,

where δ = 1/γ; note that 0 < δ < 1.

THE GENERAL CASE

Recall that b(s) = f(s) − s, where f is the pgf given by f(s) = a + _{∑i≥1} b_{i}s^{i+1} , and D = −b′(1−). Recall also that, in the supercritical case, the absorption probabilities have generating function

∑i≥1 α_{i}s^{i} = s/(1 − s) − Ds/b(s)

(α_{i} = 1 for all i ∈ C in the subcritical case). We will use the fact that b(s) = 0 has a unique solution σ on [0, 1], and that σ = 1 or 0 < σ < 1 according as D ≥ 0 or D < 0.

Lemma 1. λ = −ρb′(σ−).

Theorem 1. In the subcritical case both types of QSD exist. The Type I QSD is given by π_{j}^{(1)} = (1 − σ)σ^{j−1} (nice!), and the Type II

QSD has pgf Π^{(2)} (s) = X(σs)/X(σ−), where

X(s) = \frac{s}{b(s)} exp (−λB(s)), s < σ,

and, for s < σ, B(s) = ρ^{−1} ∫_{0}^{s} dy/b(y).

THE GENERAL CASE

The supercritical case is more delicate. We require two extra conditions:

Condition (A). The catastrophe-size distribution has finite second moment, that is, f′′(1−) < ∞ (equivalently b′′(1−) < ∞).

Condition (B). The function b can be written

b(s) = D(1 − s) + (1 − s)^{2}L((1 − s)^{−1}),

where L is slowly varying, that is, L(xt) ∼ L(x) for large t.

Theorem 2. In the supercritical case both types of QSD exist under (A) and (B). The Type I QSD is given by π_{j}^{(1)} = α_{j}/ ∑_{k≥1} α_{k},

and the Type II QSD has pgf Π^{(2)} (s) = X(s)/X(1−), where

X(s) = \frac{s}{b(s)} exp (−λB(s)), s < 1,

and, for s < 1, B(s) = ρ^{−1} ∫_{0}^{s} dy/b(y).