LIMITING CONDITIONAL DISTRIBUTIONS FOR STOCHASTIC METAPOPULATION

MODELS

by

Phil Pollett

Department of Mathematics The University of Queensland

CONDITIONAL STATE DISTRIBUTION

X(t) - state of the population at time t

We suppose that (X(t), t ≥ 0) is a continuous time Markov chain with a discrete state space S = {0} ∪ C, where 0 is the state corresponding to extinction and C comprises the remaining states.

p_{x}(t) = Pr(X(t) = x) - state probabilities

We observe the population at an arbitrary time u and extinction has not yet occurred. How can we incorporate this information?

The conditional state distribution

Evaluate the state probabilities at time u conditioned on non-extinction:

m_{x}(u) = Pr(X(u) = x|X(u) ≠ 0)

= \frac{p_{x}(u)}{1 − p_{0}(u)} , x ∈ C.

A STOCHASTIC MODEL FOR METAPOPULATIONS

There are n geographical regions (patches): N = {1, 2, …, n} Let X = (X_{1}, X_{2}, …, X_{n}), where X_{i}(t) is 1 or 0 according as patch i is occupied or not at time t (t = 0, 1, 2, …). Note that the state space is S = {0, 1}^{n}.

Q = (q_{ij}, i, j ∈ N) - Interaction matrix:

q_{ij}, for j ≠ i, is the probability that patch j will not be colonized by migration from patch i, and q_{ii} is the probability that (in the absence of immigration) patch i will become extinct.

Assume (Gyllenberg and Silvestrov) that

q_{ij} = exp (−e^{−ad_{ij}}A_{i}), i, j ∈ N ,

where d_{ij} is the distance between patches i and j (d_{ii} = 0 and d_{ij} = d_{ji}), A_{i} is the area of patch i and a(≥ 0) measures how badly individuals are at migrating.

THE TRANSITION PROBABILITIES

Assume that the various colonization processes and local extinction processes are independent. Define q_{i}(x), where x = (x_{1}, x_{2}, …, x_{n}), by

q_{j}(x) = n ; ∏ ; i=1 q_{ij}^{x_{i}}, j ∈ N , x ∈ S,

to be the probability that patch j will become extinct at the next epoch given a present configuration x.

The transition matrix P = (p(x, y), x, y ∈ S):

p(x, y) = n ; ∏ ; i=1 q_{i}(x)^{1−y_{i}}(1 − q_{i}(x))^{y_{i}}, x, y ∈ S.

Note that, since q_{i}(0) = 1, i ∈ N, state 0 = (0, 0, …, 0) (corresponding to the extinction of all patches) is an absorbing state for the chain:

p(0, y) = \begin{cases}1, & if y = 0, \\ 0, & otherwise . \end{cases}

A 5-PATCH METAPOPULATION

Patches 2, 3, 4 & 5 are equally spaced (a distance 0.1 apart), while Patch 1 is 10 times that distance away from the others (d_{1j} = 1). All patches have the same area.

PERSISTENCE OF METAPOPULATIONS: QUASI-STATIONARITY

Simulation of the 5-patch model with a = 7. The number of occupied patches is plotted against time (up to extinction at t = 728).

Do the conditional state probabilities account for the observed behaviour?

Comparison of the observed frequencies with the conditional state distribution m_{x}(t) at t = 1, 2, 5, 10. The brown bar is the proportion of time for which i patches were occupied (i = 1, 2, …, 5) during the period of the simulation. The blue bar is the distribution of the number of occupied patches evaluated using m_{x}(t).

LIMITING CONDITIONAL DISTRIBUTIONS

The observed trend has a simple theoretical explanation. Since C is a finite set, the limit

lim_{t→∞} m_{x}(t) = m_{x}

exists and defines a proper distribution m = (m_{x}, x ∈ C), called a limiting conditional distribution, and m is the left eigenvector of P_{C} (P restricted to C) corresponding to the eigenvalue, ρ_{1}, with maximal modulus (Darroch and Seneta (1965)). Note that the expected time till absorption, τ, is approximately ρ_{1}/(1 − ρ_{1}).

We can be precise about the rate of convergence by examining the eigenvalue, ρ_{2}, of P_{C} with second-largest modulus. It might not be real, and it has multiplicity κ ≥ 1 (for simplicity, suppose κ = 1). It can be shown that

m_{x}(t) = m_{x} + O(β^{t}) as t → ∞,

where β = |ρ_{2}|/ρ_{1}(< 1).

THE CONVERGENCE OF m_{x}(t) TO m_{x}

Comparison between the conditional state distribution (blue) and the limiting conditional distribution (brown)

of the number of occupied patches for the 5-patch metapopulation model with a = 7.

For the 5-patch metapopulation model with a = 7, we find that ρ_{1} ≃ 0.9979, ρ_{2} ≃ 0.6312 (real with multiplicity 1), β(= ρ_{2}/ρ_{1}) ≃ 0.6325 and τ ≃ 488.

PSEUDO-STATIONARY DISTRIBUTIONS (The method of Gyllenberg and Silvestrov)

Rationale: If we had assumed that Patch 1 (say) had a zero local extinction probability (q_{11} = 0), that patch would behave as a mainland. State 0 would no longer be accessible from all states: C would decompose into two classes, C_{0} and C_{1}, consisting of those states in C which have x_{1} = 0 and x_{1} = 1 respectively; either the process would start in C_{1} (mainland inhabited) and remain there, or, start in C_{0} (mainland uninhabited) and eventually enter either C_{1} or the absorbing state.

We identify a “quasi-mainland”, namely a single patch i with q_{ii} small (say Patch 1), and consider a sequence of processes indexed by ϵ = q_{ii}, treating ϵ as a perturbation.

PERTURBATION THEORY

Let ϵ ∈ (0, 1] (now arbitrary) and suppose that our interaction matrix depends on ϵ in the following way:

q_{ij}^{(ϵ)} = q_{ij} + ϵq̂_{ij} + ∘(ϵ) as ϵ → 0,

where

q_{ij} = lim_{ϵ→0} q_{ij}^{(ϵ)} and q̂_{ij} = lim_{ϵ→0} \frac{1}{ϵ} (q_{ij}^{(ϵ)} − q_{ij}) ,

the latter assumed to be non-negative and finite, and, that Q = (q_{ij}, i, j ∈ N) satisfies q_{11} = 0.

Then, in an obvious notation,

p^{(ϵ)}(x, y) = p(x, y) + ϵp̂(x, y) + ∘(ϵ), x, y ∈ S,

where P^{(ϵ)} = (p^{(ϵ)}(x, y), x, y ∈ S) is the transition matrix corresponding to Q^{(ϵ)} and P = (p(x, y), x, y ∈ S) is the transition matrix corresponding to Q.

THE G&S LIMITING REGIME

Let ϵ → 0 and t(= t_{ϵ}) → ∞ in such a way that ϵ t_{ϵ} → s, where 0 ≤ s ≤ ∞. Since the expected lifetime of the quasi-mainland is of order 1/ϵ, one is able to study the process on different time scales (s = 0, s = ∞, and 0 < s < ∞).

G&S showed that the limit

lim_{ϵ→0} Pr(X(t_{ϵ}) = y|X(0) = x), x, y ∈ C,

exists and is given by a mixture of the limiting probabilities π(x, y) for the (ergodic) chain generated by Q and the degenerate distribution δ(y, 0) which assigns all its mass to state 0, the mixing probability being e^{−λs}, where λ is a positive constant which is specified in terms of p̂(x, y).

A COMPARISON USING THE 5-PATCH MODEL

Comparison between the limiting conditional distribution (blue), the simulated proportions (green)

and the pseudo-stationary distribution (brown).

The disparity is marked: for this example, the two ways of analysing the model lead to quite different predictions.

THE EFFECT OF VARYING s

The effect of varying s on the pseudo-stationary distribution (blue). The brown bar represents the

simulated proportions of occupied patches.

Observe that the disparity becomes worse as the time-scale parameter s increases.

RECONCILIATION

Denote the state probabilities corresponding to P^{(ϵ)} by p^{(ϵ)}(t) = (p_{x}^{(ϵ)} (t), x ∈ S), and denote the corresponding conditonal probabilities by m_{x}^{(ϵ)} (t). Gosselin (1997) proved that

lim_{ϵ→0} lim_{t→∞} m_{x}^{(ϵ)} (t) = \begin{cases}π(x), & if x ∈ C_{1}, \\ 0, & if x ∈ C_{0}, \end{cases}

which he compared with Theorem 6.2 of G&S:

lim_{ϵ→0} m_{x}^{(ϵ)} (t_{ϵ}) = \begin{cases}π(x), & if x ∈ C_{1}, \\ 0, & if x ∈ C_{0}. \end{cases}

Thus, in the important case s = 0, the limiting conditonal distribution and the pseudostationary agree when ϵ is small.

The problem with the 5-patch model is that ϵ(= q_{11} ) ≃ 0.3679 – not sufficiently small.

REMARKS

Quasi-stationarity is a property of the model and not the means of analysing it. The 5-patch model exhibits quasi-stationarity, demonstrated emphatically using simulation, yet q_{11} is not small. The pseudo-stationary distribution does not capture this behaviour. On the other hand, the conditional state distribution m(t) does: afterall, it is the most information our model can provide at any time t given that we know extinction has not occured by time t. In cases when the convergence of m(t) to the limiting conditional distribution m is rapid, this distribution can be used instead.