Insect Outbreak

Fig. 1 Spruce Budworm are native to northern spruce and fir forests of the Eastern United States and Canada. For more information check out https://en.wikipedia.org/wiki/Spruce_budworm

Model of outbreak of spruce Budworm

Model for the budworm population dynamics is

$N'= rN \bigg (1-\frac{N}{K} \bigg) - p(N) \tag{1}$

Where, $r$ = growth rate, and $K$ =carrying capacity.
$N(t)$ = budworm population at time t. It is assumed to grow logistically.
$p(N)$= death rate due to predation (by birds etc.)
$\displaystyle p(N)= \frac{BN^2}{A^2+N^2}$ ; where constants A and B > 0.
$B$ is the product of bird population and predator efficiency (how good the predators are at killing the budworm).
$A$ is the switching value, the population at which the predators starts showing greater interest in harvesting the budworm.

In [1]:
# import all the packages required 
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
In [2]:
# plot of p(N) for a particular value of A and B, say A=10 and B=5

def pred(N,A,B):
    return np.array([(B*N**2)/(A**2+N**2)])

N = np.linspace(0.0, 80, 100)
plt.plot(N,pred(N,10.0,5.0)[0],'b-',lw=2);
plt.xlabel('$N$',fontsize=20)
plt.ylabel('$p(N)$',fontsize=20)
plt.axhline(5.0,color='k',ls='--')
plt.axvline(10.0,color='k',ls='--')
plt.ylim(0,6.0)

plt.text(-7.1,4.9,'B',fontsize=18)
plt.text(9.4,-0.8,'A',fontsize=18)
plt.show;

Fig 2. When budworms are scare there is not much predation, as $N=A$ predation continue to increase and finally saturates to to the value $B$. Birds are eating the budworm as fast as they can. In the plot above $A=10$ and $B=5.0$

Full model after incorporating $p(N)$ will be:

$\displaystyle N'= rN \bigg (1-\frac{N}{K} \bigg) - \frac{BN^2}{A^2+N^2} \tag {2}$

Nondimesionalization of the model

Let $u=N/A$, remember $N$ and $A$ has the same dimension. Using $\displaystyle\frac{du}{dt}=\frac{1 }{A}\frac{dN}{dt}$ and dividing both sides of eq (2) by $B$. We will get

$\displaystyle\frac{A}{B}\frac{du}{d t}= \frac{r}{B}Au\bigg (1-\frac{Au}{K} \bigg) - \frac{u^2}{1+u^2} $

Rescaling $t$ and introducing dimensionless quantities $\gamma$ and $k$ as follows:

$\displaystyle\tau=\frac{Bt}{A}$, $\displaystyle\gamma=\frac{rA}{B}$ and $\displaystyle k=\frac{K}{A}$

Then our model will be

$\displaystyle\frac{du}{d\tau}=\gamma u \bigg (1-\frac{u}{k} \bigg)- \frac{u^2}{1+u^2} \tag {3}$

Stability Analysis

One of the equilibria is $u^*=0$.

The other equilibria are given by the solutions of

$\displaystyle\gamma\bigg(1-\frac{u}{k} \bigg) = \bigg(\frac{u}{1+u^2}\bigg) \tag{4}$

In [3]:
# plot the left hand and right hand side of equation 4

def lhs_2(u):
    return np.array([u/(1.0+u**2.0)])

def rhs_2(u,g,k):
    return np.array([g*(1.0-u/k)])

u = np.linspace(0.0, 15, 10000)
plt.plot(u,lhs_2(u)[0],'b-',lw=2)

plt.plot(u,rhs_2(u,0.1,15.0)[0],'r--',lw=2)
intpt = np.argwhere(np.isclose(rhs_2(u,0.1,15.0)[0], lhs_2(u)[0], atol=0.0003)).reshape(-1)
plt.plot(u[intpt],lhs_2(u)[0][intpt],'go',ms=10)



plt.plot(u,rhs_2(u,0.2,15.0)[0],'r--',lw=2)
intpt = np.argwhere(np.isclose(rhs_2(u,0.2,15.0)[0], lhs_2(u)[0], atol=0.0003)).reshape(-1)
plt.plot(u[intpt],lhs_2(u)[0][intpt],'go',ms=10)



plt.plot(u,rhs_2(u,0.4,15.0)[0],'r-',lw=2)

intpt = np.argwhere(np.isclose(rhs_2(u,0.4,15.0)[0], lhs_2(u)[0], atol=0.0003)).reshape(-1)
plt.plot(u[intpt],lhs_2(u)[0][intpt],'go',ms=10)

plt.text(0.6,0.4,'a',fontsize=15)
plt.text(3.1,0.33,'b',fontsize=15)
plt.text(12.1,0.11,'c',fontsize=15)

plt.ylim(0.0,0.55)
plt.xlabel('$u$',fontsize=20);

Fig.2: In the figure above, blue line represents the right hand side of equation (4) and red lines the left hand side of equation (4). The value of $\gamma$ is point on y-axis where red lines intersect it and value of $k$ is lines where red lines intersect the x-axis. Intersection points (green dots) are the equilibria for given $\gamma$ and $k$. If we keep $k$ fixed and decrease $\gamma$ the number of intersection points (green dots/equilibria) decreases. That is number of equilibrium points decreases. The intercept of red line on the the x-axis is the value of $k$ and the the intercept of red line on the y-axis is the value of $r$

In [4]:
# Plot the equilibrium points on the phase line 
g=0.4;k=15.0
def f(u):
    return [g*u*(1-u/k)-(u**2)/(1+u**2)]

t = np.linspace(0.0, 15.1, 200)
ux=np.linspace(0.0, 15.1, 200)

plt.figure(figsize=(14,4))
plt.plot(ux,f(ux)[0],'b-',lw=4)
plt.axhline(0.0,color='k')
plt.xlim(-0.2,13.5);plt.ylim(-0.2,0.6)
intpt = np.argwhere(np.isclose(f(ux)[0],f(ux)[0]*0.0, atol=0.01)).reshape(-1)
plt.plot(ux[intpt],f(ux)[0][intpt],'go',ms=9)
plt.plot(0.0,0.0,'ro',ms=9)
plt.xlabel('$u$',fontsize=24);
plt.ylabel(r'$\frac{du}{d\tau}$',fontsize=24,rotation='horizontal');

plt.text(0.5,0.03,'a',fontsize=15)
plt.text(2.4,0.03,'b',fontsize=15)
plt.text(11.9,0.03,'c',fontsize=15)

plt.show()

Fig.3 red dot is the equilibrium corresponding to $u^*=0$ and green dots are the equilibrium corresponding to the Fig. 2

In [5]:
# Plot the equilibrium points on the phase line with arrows and 
#with filled/empty circle for stable/unstable equilibrium 

plt.figure(figsize=(14,4))
plt.plot(ux,f(ux)[0],'b-',lw=2)
plt.axhline(0.0,color='k')
plt.xlim(-0.2,13.5);plt.ylim(-0.2,0.6)
intpt = np.argwhere(np.isclose(f(ux)[0],f(ux)[0]*0.0, atol=0.01)).reshape(-1)

plt.plot(ux[intpt[1]],f(ux)[0][intpt[1]],'ko',ms=12)
plt.quiver(1.5, 0.0, -1, 0, pivot='mid')

plt.plot(ux[intpt[2]],f(ux)[0][intpt[2]],'o',color='none',markeredgewidth=1.5,ms=10)
plt.quiver(4.5, 0.0, 1, 0, pivot='mid')

plt.plot(ux[intpt[3]],f(ux)[0][intpt[3]],'ko',ms=12)
plt.quiver(12.6, 0.0, -1, 0, pivot='mid')

plt.plot(0.0,0.0,'o',color='none',markeredgewidth=1.5,ms=12)
plt.quiver(-0.07, 0.0, 1, 0, pivot='mid',width=0.003,)

plt.xlabel('$u$',fontsize=24);
plt.ylabel(r'$\frac{du}{d\tau}$',fontsize=24,rotation='horizontal');

plt.text(0.5,0.03,'a',fontsize=15)
plt.text(2.4,0.03,'b',fontsize=15)
plt.text(11.9,0.03,'c',fontsize=15)


# save the values of stable point a and stable point c for later use 
a_point=np.array([ux[intpt[1]],f(ux)[0][intpt[1]]])
b_point=np.array([ux[intpt[2]],f(ux)[0][intpt[2]]])
c_point=np.array([ux[intpt[3]],f(ux)[0][intpt[3]]])
                  

    
plt.show()

Fig. 4: Equilibrium $a$ and $c$ are stable, whereas $b$ is unstable. The equilibrium corresponding to $u^{*}=0$ is unstable.

Point a is called the refuge level. Point b is called the out-break level. Pest control would like to keep the population at stable equilibrium a away from point c. Whether final dynamics will end up at a or c will depend on the initial conditions as we will see below.

In [6]:
# Plot trajectories of the dynamics for different initial conditions of u
from scipy.integrate import odeint


def fun(u, t, g,k): 
    f = g*u*(1-u/k)-(u**2)/(1+u**2)
    return f


plt.figure(figsize=(10,6))

# Times at which the solution is to be computed.
t = np.linspace(0, 130, 500)

set_of_initial_conditions=np.array([0.001,0.01,0.1,0.2,0.6,1.0,2.0,2.6,3.0,5.0,6.0,7.0,10,11,13,15,20])

for uinit in set_of_initial_conditions:
    # Initial condition
    u0 =uinit

    # Parameter value to use in `fun`.
    g=0.4;k=15.0;

    # Solve the equation.
    u = odeint(fun, u0, t, args=(g,k,))

    
    plt.plot(t, u[:,0],lw=2)

    
    
plt.plot(t[-1],a_point[0],'ko',ms=15)
    
plt.plot(t[-1],c_point[0],'ko',ms=15)
plt.plot(t[-1],b_point[0],'o',color='none',markeredgewidth=1.5,ms=12)    
 

plt.xlabel('t',fontsize=25)
plt.ylabel('u',fontsize=25)
plt.xlim(0.0,140);plt.ylim(0.0,15.5)

plt.text(133,0.2,'a',fontsize=22)
plt.text(133,2.35,'b',fontsize=22)
plt.text(133,11.5,'c',fontsize=22)


plt.show()

Fig 5: Different trajectories that system can take when we start with different initial conditions. Observe that all the trajectory finally merge into either stable equilibrium a or c.

Fig 5 shows the importance of the intial condition in determining whether the end state will a or c.