Notebook

2.5.1 Summary¶

This notebook provides an introduction to nonlinear dynamics using a well-known model for the preditor-prey interaction of Snowshoe Hare and Canadian Lynx. Topics include limit cycles, the existence of multiple steady states, and simple phase plane analysis using nullclines. This notebook can be displayed as a slide presentation.

2.5.2 Introduction¶

Snowshoe hare (Lepus americanus) are the primary food for the Canadian lynx (Lynx canadensis) in the Northern boreal forests of North America. When hare are abundant, Lynx will eat hare about two every three days almost to the complete exclusion of other foods. As a consequence, the population dynamics of the two species are closely linked.

Canadian Lynx	Snowshoe Hare

kdee64 (Keith Williams) CC BY 2.0, via Wikimedia Commons	D. Gordon E. Robertson CC BY-SA 3.0, via Wikimedia Commons

It has been known for over a century that the populations of the two species vary dramatically in cycles of 8 to 11 year duration. This chart, for example, shows pelt-trading data taken from the Hudson's Bay Company (from MacLulich, 1937. See important notes on this data in Stenseth, 1997)

(CNX OpenStax CC BY 4.0, via Wikimedia Commons)

The actual cause of the cycling is still a matter of scientific inquiry. Hypotheses include the inherent instability of the preditor-prey dynamics, the dynamics of a more complex food web, and the role of climate (see Zhang, 2007). The discussion in this notebook addresses the preditor-prey dynamics.

2.5.3 Historical Data¶

A digitized version of the historical data is available from D. R. Hundley at Whitman College. The following cell reads the data from the url, imports it into a pandas dataframe, and creates a plot.

In [20]:

%matplotlib inline
import pandas as pd

url = 'http://people.whitman.edu/~hundledr/courses/M250F03/LynxHare.txt'
df = pd.read_csv(url, delim_whitespace=True, header=None, index_col=0)
df.index.name = 'Year'
df.columns = ['Hare', 'Lynx']
df.plot(figsize=(10,6), grid=True)

Out[20]:

<matplotlib.axes._subplots.AxesSubplot at 0x11dfab210>

2.5.4 Population Dynamics¶

2.5.4.1 Model Equations¶

The model equatons describe the time rate of change of the population densities of hare ($H$) and lynx ($L$). Each is the difference between the birth and death rate. The death rate of hare is coupled to the population density of lynx. The birth rate of lynx is a simple multiple of the death rate of hare.

$$\begin{align*}\frac{dH}{dt} & = \underbrace{rH\left(1-\frac{H}{k}\right)}_{Hare Birth Rate}-\underbrace{\frac{aHL}{c+H}}_{Hare Death Rate}\\ \frac{dL}{dt} & = \underbrace{a\frac{bHL}{c+H}}_{Lynx Birth Rate}-\underbrace{dL}_{Lynx Death Rate} \end{align*}$$

2.5.4.2 Parameter Values¶

Parameter	Symbol	Value
Lynx/Hare Predation Rate	$a$	3.2
Lynx/Hare Conversion	$b$	0.6
Lynx/Hare Michaelis Constant	$c$	50
Lynx Death Rate	$d$	0.56
Hare Carrying Capacity	$k$	125
Hare Reproduction Rate	$r$	1.6

2.5.5 Simulation using the scipy.odeint()¶

2.5.5.1 Step 1. Initialization¶

The SciPy library includes functions for integrating differential equations. Of these, the function odeint provides an easy-to-use general purpose algorithm well suited to this type of problem.

2.5.5.2 Step 2. Establish Parameter Values¶

Set global default values for the parameters

2.5.5.3 Step 3. Write function for the RHS of the Differential Equations¶

deriv is a function that returns a two element list containting values for the derivatives of $H$ and $L$. The first argument is a two element list with values of $H$ and $L$, followed by the current time $t$.

$$\begin{align*} \frac{dH}{dt} & = r H \left(1-\frac{H}{k}\right) - \frac{a H L}{c + H} \\ \frac{dL}{dt} & = b\frac{a H L}{c + H} - dL \end{align*}$$

2.5.5.4 Step 4. Choose Time Grid, Initial Conditions, and Integrate¶

In [25]:

# perform simulation
t = np.linspace(0,70,500)       # time grid
IC = [20,20]                    # initial conditions for H and L
sol = odeint(deriv,IC,t)        # compute solution
H,L = sol.transpose()           # unpack solution 

2.5.5.5 Step 5. Visualize and Analyze the Solution¶

For this choice of parameters and initial conditions, the Hare/Lynx population exhibits sustained oscillations.

2.5.5.5.1 Phase Plane¶

In [27]:

plt.figure(figsize=(13,4))
plt.subplot(1,2,1)
plt.plot(t,H,t,L)
plt.title('Hare/Lynx Population Dynamics')
plt.xlabel('Year')
plt.legend(['Hare','Lynx'])
plt.grid(True)

plt.subplot(1,2,2)
plt.plot(H,L)
plt.title('Hare/Lynx Phase Plot')
plt.ylabel('Lynx')
plt.xlabel('Hare')
plt.grid(True)

2.5.6 Nullclines¶

Nullclines are the points in the phase plane where the derivatives are equal to zero.

The nullclines for hare are where

$$\frac{dH}{dt} = 0 \implies \begin{cases} \begin{align*} H^* & = 0 \\ \\ L^* & = \frac{r}{a}\left(c+H\right)\left(1-\frac{H}{k}\right) \end{align*} \end{cases}$$

The nullclines for Lynx are where

$$\frac{dL}{dt} = 0 \implies \begin{cases} \begin{align*} L^* & = 0 \\ \\ H^* & = \frac{c d}{a b - d} \end{align*} \end{cases}$$

For convenience, we create a function to plots the nullclines and steady states that occur where the nullclines intersect.

In [28]:

def plotNullclines():
    # nullcline dH/dt = 0
    Hp = np.linspace(0,k)
    Lp = (r/a)*(c+Hp)*(1-Hp/k)
    plt.plot(Hp,Lp,'b')
    plt.ylim(0,130)
    plt.xlim(0,150)
    
    # nullcline dL/dt = 0
    Hd = c*d/(a*b-d)
    plt.plot([Hd,Hd],plt.ylim(),'r')

    # additional nullclines
    plt.plot([0,0],plt.ylim(),'b')
    plt.plot(plt.xlim(),[0,0],'r')

    # steady states
    Hss = c*d/(a*b-d)
    Lss = r*(1-Hss/k)*(c+Hss)/a
    plt.plot([0,k,Hss],[0,0,Lss],'r.',ms=20)

    plt.xlabel('Hare')
    plt.ylabel('Lynx')
    plt.legend(['dH/dt = 0','dL/dt = 0'])  

Here's a plot of the nullclines for the default parameter values. The steady states correspond to

No Hare, and no Lynx.
Hare population at the carry capacity of the environment, and no Lynx
Coexistence of Hare and Lynx.

Visualization of the nullclines give us some insight into how the Hare and Lynx populations depend on the model parameters. Here we look at how the nullclines depend on the Hare/Lynx predation rate $a$.

In [30]:

from ipywidgets import interact

def sim(aslider= 3.2):
    global a
    a = aslider
    plt.xlim(0,150)
    plt.ylim(0,130)
    plotNullclines()
    
interact(sim,aslider=(1.25,4,.01))

Out[30]:

<function __main__.sim(aslider=3.2)>

2.5.7 Interactive Simulation¶

2.5.7.1 Visualization Function¶

The visualization function for this example accepts a list of time values, values of $H$ and $L$, and model parameters. The model parameters are needed to plot nullclines and steady states on the phase plane.

In [31]:

# visualization
def HLPlot(t,H,L):
    # time axis
    plt.figure(figsize=(12,5))
    plt.subplot(1,2,1)
    plt.plot(t,H)
    plt.plot(t,L)
    plt.xlabel('Time [years]')
    plt.ylabel('Population Density')
    plt.legend(['Hare','Lynx'],loc='upper left')
    
    # phase plane
    plt.subplot(1,2,2)
    plt.plot(H,L)
    plt.xlim(0,150)
    plt.ylim(0,130)
    plotNullclines()

2.5.7.2 Simulation Function¶

An additional function is created to encapsulate the entire process of solving the model and displaying the solution. The function takes arguments specifing the initial values of $H$ and $L$, and a value of the parameter $a$. These argument

In [32]:

# default parameter values
a = 3.2
b = 0.6
c = 50
d = 0.56
k = 125
r = 1.6

# perform simulation
t = np.linspace(0,70,500)

def LynxHare(H=20,L=20,aslider=3.2):
    IC = [H,L]
    global a
    a = aslider
    sol = odeint(deriv,IC,t)
    HLPlot(t,sol[:,0],sol[:,1])

Use the aslider to adjust values of the Hare/Lynx interaction. Can you indentify stable and unstable steady states?

2.5.8 Stability of a Steady State¶

2.5.8.1 1. Unstable Focus¶

Any displacement from an unstable focus leads to a trajectory that spirals away from the steady state.

2.5.8.2 2. Stable Focus¶

Small displacements from a stable focus results in trajectories that spiral back towards the steady state.

2.5.8.3 3. Stable and Unstable Nodes¶

Displacements from a steady state either move towards (stable) or away from (unstable) nodes without the spiral structure of a focus.

2.5.9 Summary¶

Hope you enjoyed this brief introduction to the modeling of a small food web. This is a fascinating field with many important and unanswered questions. Recent examples in the research literature are here and here.

What you should learn from this notebook:

How to simulate small systems of nonlinear ODEs.
How to plot trajectories in a phase plane.
How to plot the nullclines of two differential equations with constant parameters.
Solve systems for multiple steady states.
Recognize limit cycles, steady-states, stable and unstable foci, stable and unstable nodes.

2.5.10 Suggested Exercise¶

Explore the impact of the parameter $a$ on the nature of the solution. $a$ is proporational to the success of the Lynx hunting the Hare. What happens when the value is low? high? Can you see the transitions from conditions when the Lynx done't survive, the emergence of a stable coexistence steady-state, and finally the emergence of a stable limit cycle?

2.5 Hare and Lynx Population Dynamics¶