# setup
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import scipy
import scipy.integrate
# setup: has to be in a separate cell for some reason
plt.rcParams['figure.figsize'] = [12, 7]
Suppose we are modeling the population sizes
of two interacting species,
wihch we creatively call "type 0
" and "type 1
".
The dynamics are:
Note: for simplicity, we've kept death and reproduction coupled, which might seem wierd, but imagine annual plants.
How do we expect this to behave? Let's write down the expected dynamics: E[Nt+1|Nt=N,Mt=M]=N(1−p0)+p0Nλ0exp(−(N+a0M)/K0)E[Mt+1|Nt=N,Mt=M]=M(1−p1)+p1Mλ1exp(−(M+a1N)/K1).
The amount that N changes per unit time is Nt+1−Nt. This suggests the following differential equation: dNdt=p0Nt(λ0exp(−Nt+a0MtK0)−1)dMdt=p1Mt(λ1exp(−Mt+a1NtK1)−1)
We'll be using these again, so we'll denote this by F:R2↦R2: F(N,M)=(p0N(λ0exp(−(N+a0M)/K0)−1),p1M(λ1exp(−(M+a1N)/K1)−1)).
def run_sim_2d(N0, M0, gen_fn, ngens, dtype='int', **kwargs):
"""
N0, M0: initial population size arrays for the two populations
gen_fn: function that calculates the next generation of population sizes
given the current population sizes. Takes in two arrays and other arguments
and outputs a tuple of length two of the next generation's population sizes.
ngens: number of generations to do this for
dtype: the numpy type (eg either 'int' or 'float') that gen_fn returns
**kwargs : the other parameters for gen_fn
returns: a tuple of arrays, each of size (ngens, len(N0)),
with columns corresponding to replicates and rows corresponding to time points
"""
assert(len(N0) == len(M0))
N = np.empty((ngens, len(N0)), dtype=dtype)
N[0, :] = N0
M = np.empty((ngens, len(M0)), dtype=dtype)
M[0, :] = M0
for t in range(1, ngens):
n, m = gen_fn(N[t-1, :], M[t-1, :], **kwargs)
assert(n.dtype == dtype and m.dtype == dtype)
N[t, :], M[t, :] = n, m
return N, M
Here's what this looks like:
def lv_eqn(N, M, lam, p, K, a):
next_N = N + p[0] * N * (lam[0] * np.exp(-(N + a[0] * M)/K[0]) - 1)
next_M = M + p[1] * M * (lam[1] * np.exp(-(M + a[1] * N)/K[1]) - 1)
return next_N, next_M
def lv_gen(N, M, lam, p, K, a):
assert(len(M) == len(N))
N_dies = np.random.binomial(N, p[0])
M_dies = np.random.binomial(M, p[1])
N_repro = np.random.poisson(N_dies * lam[0] * np.exp(-((N + a[0]*M)/K[0])))
M_repro = np.random.poisson(M_dies * lam[1] * np.exp(-((M + a[1]*N)/K[1])))
next_N = N - N_dies + N_repro
next_M = M - M_dies + M_repro
return next_N, next_M
# we'll use **kwargs to avoid re-typing these over and over
lv_args = {
'lam' : [1.2, 1.1], # fecundity
'p' : [0.5, 0.9], # death prob
'K' : [2000, 5000], # carrying capacity
'a' : [0.1, 0.1] # encounter rate
}
# random simulation
N, M = run_sim_2d([500], [100], lv_gen, 120, dtype='int',
**lv_args)
# deterministic
tN, tM = run_sim_2d([500], [100], lv_eqn, 120, dtype='float',
**lv_args)
fig, ax = plt.subplots()
ax.plot(np.column_stack([N, M]))
ax.plot(np.column_stack([tN, tM]), linestyle="--")
Equilibria occur when F(N,M)=0.
First note that if M=0 then the model is one-dimensional, and so there are equilibria on the boundary. These occur when E[Nt+1]=N, which happens only if the average number of new offspring for a given parent is 1, i.e., F(N0,0)=0⇒N0=K0log(λ0)F(0,M0)=0⇒M0=K1log(λ1).
For an internal equilibrium, i.e., one where both co-exists, we require that λ0exp(−(N+a0M)/K0)=1, and this happens only if (N+a0M)/K0=logλ0.
Writing out the equations for both species, equilibrium occurs when N+a0M=K0log(λ0)a1N+M=K1log(λ1) which is solved by N∗=K0log(λ0)−a0K1log(λ1)1−a0a1M∗=K1log(λ1)−a1K0log(λ0)1−a0a1.
These only make sense if they are positive... but let's have a look at the pictures, first.
def lv_equil(lam, p, K, a):
bequil = np.array([[K[0] * np.log(lam[0]), 0],
[0, K[1] * np.log(lam[1])]])
equil = np.array([K[0] * np.log(lam[0]) - a[0] * K[1] * np.log(lam[1]),
K[1] * np.log(lam[1]) - a[1] * K[0] * np.log(lam[0])]) / (1 - a[0] * a[1])
return np.row_stack([equil, bequil])
def phase_plot(genfn, xlim, ylim, nx, ny, scale=1, figax=None, **kwargs):
xstep = int((xlim[1] - xlim[0]) / nx)
ystep = int((ylim[1] - ylim[0]) / ny)
X, Y = np.meshgrid(range(xlim[0], xlim[1], xstep), range(ylim[0], ylim[1], ystep))
X.shape = Y.shape = (np.prod(X.shape),)
U, V = genfn(X, Y, **kwargs)
if figax is None:
figax = plt.subplots()
fig, ax = figax
ax.quiver(X, Y, U-X, V-Y, angles='xy', scale_units='xy', scale=scale)
return fig, ax
def do_plot(lv_args, init = [[500], [100]], nsteps=120, figax=None):
equil = lv_equil(**lv_args)
tN, tM = run_sim_2d(init[0], init[1], lv_eqn, nsteps, dtype='float',
**lv_args)
fig, ax= phase_plot(lv_eqn, xlim=[0, 600], ylim=[0, 600], nx=20, ny=20,
figax=figax, **lv_args)
ax.plot(tN, tM)
ax.scatter(tN, tM)
ax.scatter(equil[0,0], equil[0,1], s=1000, c='g', alpha=0.5)
ax.scatter(equil[1:,0], equil[1:,1], s=1000, c='r', alpha=0.5)
ax.axis('equal')
ax.set_xlabel("N")
ax.set_ylabel("M")
ax.set_aspect(1.0)
return fig, ax
lv_args = {
'lam' : [1.2, 1.1], # fecundity
'p' : [0.5, 0.9], # death prob
'K' : [2000, 5000], # carrying capacity
'a' : [0.1, 0.1] # encounter rate
}
fig, ax = do_plot(lv_args)
And, here is a stochastic simulation trajectory. (Run again to add more.)
# random simulation
N, M = run_sim_2d([500], [100], lv_gen, 120, dtype='int', **lv_args)
ax.plot(N, M)
ax.scatter(N, M)
fig
A nice way to look at phase plots is to draw on them the isoclines, i.e., the lines along which each variable does not change. We saw above that N does not change if N+a0M=K0log(λ0), and M does not change if a1N+M=K1log(λ1). These are just straight lines in phase space.
def plot_isoclines(ax, lam, p, K, a, Nmax=600):
xx = np.array([0.0, Nmax])
ax.plot(K[0] * np.log(lam[0]) - xx * a[0], xx)
ax.plot(xx, K[1] * np.log(lam[1]) - xx * a[1])
plot_isoclines(ax, **lv_args)
fig
But: here's a different set of parameters:
lv_args = {
'lam' : [1.2, 1.1], # fecundity
'p' : [0.5, 0.9], # death prob
'K' : [2000, 5000], # carrying capacity
'a' : [0.75, 0.75] # encounter rate
}
plt.rcParams['figure.figsize'] = [30, 20]
fig, axes = plt.subplots(2, 3)
for a, ax in zip([0.3, 0.4, 0.5, 0.6, 0.7, 0.8], axes.flatten()):
lv_args['a'] = [a, a]
fig, ax = do_plot(lv_args, figax=(fig, ax), nsteps=1000)
plot_isoclines(ax, **lv_args)
Recall the internal equilibrium of this system is N∗=K0log(λ0)−a0K1log(λ1)1−a0a1M∗=K1log(λ1)−a1K0log(λ0)1−a0a1.
In the case that a0a1<1, these are positive iff 1a1K1log(λ1)>K0log(λ0)>a0K1log(λ1)
What happens if they are not positive? Since growth without bound is not possible, the dynamics approach a boundary. But, which one? What we need to know is whether the trajectories started a little away from the equilbria tend to get closer or farther away. For instance, if FM(N0,ϵ)>0, then even small numbers of the second species will grow when the first species is at its boundary equilibrium.
To find out, we look at the gradient of F, which we denote H: HNN=ddNFN(N,M)=p0{(1−NK0)λ0exp(−(N+a0M)/K0)−1},HNM=ddNFM(N,M)=−p0a0NK0λ0exp(−(N+a0M)/K0),HMN=ddNFM(N,M)=−p1a1MK1λ1exp(−(M+a1N)/K1),HMM=ddMFM(N,M)=p1{(1−MK1)λ1exp(−(M+a1N)/K1)−1}.
So - FM(N0,ϵ)≈ϵHMM(N0,0), and since N0=K0logλ0, HMM(N0,0)=p1{λ1exp(−a1log(λ0)K0/K1)−1}=p1{λ1λ−a1K0/K10−1}. In other words, the second species, M, can invade the equilibrium of the first, N, if K1log(λ1)>a1K0log(λ0).
def lv_H(lam, p, K, a):
def H(N, M):
HNN = p[0] * ( (1 - N/K[0]) * lam[0] * np.exp(- (N + a[0] * M)/ K[0]) - 1)
HNM = - p[0] * a[0] * N / K[0] * lam[0] * np.exp(- (N + a[0] * M)/ K[0])
HMN = - p[1] * a[1] * M / K[1] * lam[1] * np.exp(- (M + a[1] * N)/ K[1])
HMM = p[1] * ( (1 - M/K[1]) * lam[1] * np.exp(- (M + a[1] * N)/ K[1]) - 1)
return np.array([[HNN, HNM], [HMN, HMM]])
return H
H = lv_H(**lv_args)
for eq in lv_equil(**lv_args):
print("Equilibrium:", eq)
print(" with H = \n", H(*eq))