In [0]:
import numpy as np
from matplotlib import pyplot as plt
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
from scipy.integrate import odeint

mpl.rcParams['legend.framealpha'] = 0.8
mpl.rcParams['legend.facecolor'] = 'w'
mpl.rcParams['legend.frameon'] = True

Basic Lotka-Volterra, Predator-Prey

This is an exploration of non-linear differential equation modeling in Python, using the Lotka-Volterra equation as an initial example.

Variables

$x \geq 0$ is the prey population (prey)

$y \geq 0$ is the predator population (predator)

$a > 0$ is the natural reproduction rate of prey (prey_reproduction)

$b > 0$ is the factor at which predators meet and kill prey (prey_predator_victim)

$c > 0$ is the natural mortality rate of predators (predator_mortality)

$d > 0$ is the factor relating the predator reproduction rate to the abundance of prey (predator_reproduction)

Differential Equations

$\frac{dx}{dt} = ax - bxy$ describes the change in prey population over time, with $ax$ representing natural exponential population growth and $-bxy$ representing deaths to predation

$\frac{dy}{dt} = dxy - cy$ describes the change in predator population over time, with $-cy$ representing natural exponential population decay and $dxy$ representing exponential population growth dependant on the abundance of prey

In [13]:
def model_population_changes(populations, time, prey_reproduction, prey_predator_victim, predator_mortality, predator_reproduction):
    prey = populations[0]
    predator = populations[1]
    dPrey = prey_reproduction * prey - prey_predator_victim * prey * predator
    dPredator = predator_reproduction * prey * predator - predator_mortality * predator
    
    return [dPrey, dPredator]

initial_populations = [
    120,    # prey
    8       # predators
]

params = (
    1.0,    # a = prey_reproduction
    0.1,    # b = prey_predator_victim
    1.0,    # c = predator_mortality
    0.01    # d = predator_reproduction
)

times = np.linspace(0, 20, 1000)

def detail_plot(initial_populations=initial_populations, params=params, times=times, title='', show=True):
    out = odeint(func=model_population_changes,y0=initial_populations,t=times,args=params)

    prey, predator = out.T

    fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(15, 5))
    plt.suptitle('{} (x={}, y={}, a={}, b={}, c={}, d={} )'.format(title, *initial_populations, *params))

    plt.subplot(121)
    plt.title('Populations / Time')
    plt.plot(times, prey, 'c', label='Prey (x)')
    plt.plot(times, predator, 'r', label='Predator (y)')
    plt.xlabel('Time')
    plt.ylabel('Population')
    plt.legend()

    axes[1].set_title('Phase Plot')
    axes[1].plot(prey, predator, color='m', alpha=0.5, zorder=1)
    
    # Plot the initial population with an arrow showing the direction of the phase change
    if prey[1] != prey[0] and predator[1] != predator[0]:
        axes[1].quiver(initial_populations[0], initial_populations[1], prey[1] - prey[0], predator[1] - predator[0], angles='xy', zorder=2) 

    axes[1].set_xlabel('Prey (x)')
    axes[1].set_ylabel('Predator (y)')
    axes[1].set_xbound(0)
    axes[1].set_ybound(0)
    
    if show: plt.show()
        
    return axes
    
detail_plot();

Equilibrium

By setting the differential equations each to zero, I found two equilibrium solutions:

Extinction: $x = y = 0$

Dependent on parameters: $x = \frac{c}{d}, y = \frac{a}{b}$

Using these formulae I can set the populations to the requisite starting values and the populations will stay constant.

Different Time Step Sizes

I played with different type step sizes using np.linspace rather than range. It was quickly evident that with large enough time steps, the model became odd and inaccurate. Also, the oscillation was less stable with large step sizes. One interesting side effect is that I learned how to make a Spirograph out of a Phase Plot.

In [14]:
detail_plot(times=np.linspace(0, 100, 50), title='Large (101/50) Time Steps')
detail_plot(times=np.linspace(0, 100, 10000), title='Small (101/10000) Time Steps');
In [15]:
equilibrium_populations = [
    params[2] / params[3], # Prey
    params[0] / params[1]  # Predators
]
axes = detail_plot(initial_populations=equilibrium_populations, show=False)
axes[1].plot(*equilibrium_populations, 'bD', label='Equilibrium')
axes[1].legend(frameon=True, framealpha=0.7, facecolor='w')
print(equilibrium_populations)
axes[1].set_xbound(0)
axes[1].set_ybound(0)
plt.show()
[100.0, 10.0]

Parameter Exploration

I then decided to experiment with plotting different parameter values but keeping the equilibrium population as the initial population. After playing with a few different visualizations I decided to use small-multiples to show variation in two different parameters, $a$ and $b$. I stuck with phase plots, because suprisingly, although these were initially very opaque (this is the first time I've seen the concept), after staring at these all day I actually feel that I get more useful information out of these.

I did however have to use the plt.quiver method to plot an arrow, as I appeciate knowing the direction the phase plot moves with the progress of time. Because of that, I noticed that if the initial predator population is higher than that of equilibrium, the prey population starts in decline. If the initial predator population is lower, the prey population starts in incline. This changes the oscillation direction at the outset of observation, as evidenced by the left or right facing arrows for the different phase plots.

In [16]:
ncols=3
nrows=3
fig, axes = plt.subplots(nrows=nrows, ncols=ncols, sharex=True, sharey=True, figsize=(15,10))
fig.suptitle('Phase Plots with Varying Prey Reproduction and Prey Victim Parameters')

prey_reproductions = np.linspace(params[0], params[0]*2, nrows)
prey_victims = np.linspace(params[1], params[1]*2, ncols)

for i in range(nrows):
    for j in range(ncols):
        params = (prey_reproductions[i], prey_victims[j], params[2], params[3])
        out = odeint(func=model_population_changes, y0=equilibrium_populations, t=times, args=params, mxstep=100000)
        prey, predator = out.T
        
        # Plot the analytical solution of dx/dt = 0, dy/dt = 0 (x=c/d, y=a/b)
        # This is an equilibrium point
        axes[i, j].plot((params[2] / params[3]), (params[0] / params[1]), 'bD', label='Equilibrium')
          
        axes[i, j].set_title('a={}, b={}'.format(round(params[0], 5), round(params[1], 5)))
        axes[i, j].set_xlabel('Prey')
        axes[i, j].set_ylabel('Predators')
        axes[i, j].label_outer()
        axes[i, j].plot(prey, predator, alpha=0.5)
        
        # Plot the initial population with an arrow showing the direction of the phase change
        if prey[1] != prey[0] and predator[1] != predator[0]:
            axes[i, j].quiver(equilibrium_populations[0], equilibrium_populations[1], prey[1] - prey[0], predator[1] - predator[0], angles='xy')        

        axes[i, j].set_ybound(lower=0)
        axes[i, j].set_xbound(lower=0)
        axes[i, j].grid()

plt.show()

Advanced Predator-Prey with Super-Predator

To model an additional layer of complexity I added a super-predator that consumes both normal predators and prey. I considered what the effect would be on each existing differential equation as well as what their own differential equation might be. Clearly, the super-predator population would have a negative impact on the prey and predator populations, just as the predator population has a negative impact on prey. The super-predator population would also naturally drop without the existance of food, like the predator's. I also decided to model the super-predator reproduction rate as independent factors of the abundance of prey and the abundance of predators, because the super-predators could reproduce even if only one food source existed (predator or prey), and super-predators may get more "nutrients" from predators than prey or vice versa.

Variables

$x \geq 0$ is the prey population (prey)

$y \geq 0$ is the predator population (predator)

$z \geq 0$ is the super-predator population (spredator)

$a > 0$ is the natural reproduction rate of prey (prey_reproduction)

$b > 0$ is the factor at which predators meet and kill prey (prey_predator_victim)

$c > 0$ is the natural mortality rate of predators (predator_mortality)

$d > 0$ is the factor relating the predator reproduction rate to the abundance of prey (predator_reproduction)

$e > 0$ is the factor at which super-predators meet and kill prey (prey_spredator_victim)

$f > 0$ is the factor at which super-predators meet and kill predators (predator_spredator_victim)

$g > 0$ is the natural mortality rate of super-predators (spredator_mortality)

$h > 0$ is the factor relating the super-predator reproduction rate to the abundance of prey (spredator_prey_reproduction)

$i > 0$ is the factor relating the super-predator reproduction rate to the abundance of predators (spredator_predator_reproduction)

Differential Equations

$\frac{dx}{dt} = ax - bxy - exz$ describes the change in prey population over time, with $ax$ representing natural exponential population growth, $-bxy$ representing deaths to predation from predators, and $-exz$ representing deaths to predation from super-predators

$\frac{dy}{dt} = dxy - cy - fyz$ describes the change in predator population over time, with $-cy$ representing natural exponential population decay, $dxy$ representing exponential population growth dependant on the abundance of prey, and $-fyz$ representing deaths to predation from super-predators

$\frac{dz}{dt} = hxz + iyz - gz$ describes the change in super-predator population over time, with $hxz$ representing exponential population growth dependant on the abundance of prey, $iyz$ representing exponential population growth dependant on the abundance of predators, and $-gz$ representing natural exponential population decay.

In [0]:
from matplotlib import animation, rc
from IPython.display import HTML
rc('animation', html='html5')

def model_population_changes_super(populations,
                                   time,
                                   prey_reproduction,
                                   prey_predator_victim,
                                   predator_mortality,
                                   predator_reproduction,
                                   prey_spredator_victim,
                                   predator_spredator_victim,
                                   spredator_mortality,
                                   spredator_prey_reproduction,
                                   spredator_predator_reproduction):
    prey = populations[0]
    predator = populations[1]
    spredator = populations[2]
    
    dPrey = prey_reproduction * prey - prey_predator_victim * prey * predator - prey_spredator_victim * prey * spredator
    dPredator = predator_reproduction * prey * predator - predator_mortality * predator - predator_spredator_victim * predator * spredator
    dSPredator = spredator_prey_reproduction * prey * spredator + spredator_predator_reproduction * predator * spredator - spredator_mortality * spredator

    return [dPrey, dPredator, dSPredator]

initial_populations_super = [
    100,    # prey
    50,     # predators
    0,      # super-predators
]

params_super = (
    1.0,    # a = prey_reproduction
    0.02,   # b = prey_predator_victim
    1.0,    # c = predator_mortality
    0.01,   # d = predator_reproduction
    0.01,   # e = prey_spredator_victim
    0.01,   # f = predator_spredator_victim
    0.0,    # g = spredator_mortality
    0.0,    # h = spredator_prey_reproduction
    0.0     # i = spredator_predator_reproduction
)

times_super = np.linspace(0, 50, 5000)

def detail_plot_super(initial_populations=initial_populations_super, times=times_super, params=params_super):
    out = odeint(func=model_population_changes_super,y0=initial_populations,t=times,args=params)

    prey, predator, spredator = out.T

    plt.figure(figsize=(15, 5))
    axes = plt.subplot(121)
    plt.title('Populations / Time')
    plt.plot(times, prey, 'c', label='Prey (x)')
    plt.plot(times, predator, 'r', label='Predator (y)')
    plt.plot(times, spredator, 'g', label='Super Predator (z)')
    plt.xlabel('Time')
    plt.ylabel('Population')
    axes.set_ybound(0)
    plt.legend()

    axes = plt.subplot(122, projection='3d')
    plt.plot(prey, predator, spredator, 'c')
    if min(prey) == max(prey) and min(predator) == max(predator) and min(spredator) == max(spredator):
        plt.plot(prey, predator, spredator, 'cD', label='Equilibrium')
        plt.legend()
    axes.set_xlabel('Prey')
    axes.set_ylabel('Predators')
    axes.set_zlabel('Super Predators')
    plt.title('Phase Plot')

    plt.show()

It became clear that mathematical analysis of the differential equations to get an equilibrium result is exponentially more difficult as you add more actors. Rather than the directly balancing negative feedback loops as seen in the 2-agent model, here there are competing predators and it becomes necessary to find a specific equilibrium point for the system given an initial population or a set of parameters. This equilibrium could either be static or oscillating. This differs from the 2-agent model where all initial conditions (given enough model fidelity) result in at least an oscillating system.

I did find some interesting information though by trying different variations of parameters and initial populations. I started my exploration with the assumption that the super predator model with an initial super predator population of $z = 0$ should be no different than the normal Lotka-Volterra model presented before, which is an assumption supported by my 3-agent model equations. By starting with initial equlibrium conditions of the previous model (this time using parameters supporting equilibrium populations of $x = 100, y = 50$) I was able to experiment with various super predator parameters and see how they effected the model.

Equilibrium State - No Super Predators

In [18]:
detail_plot_super()

Stable Oscillation - Super Predators Unaffected by Victim Populations

If super predators are unaffected by the prey and predator populations and don't die or reproduce themselves, then they create a simple stable oscillation pattern. This is shown by a 3-dimensional phase plot where the phase inhabits a plane where the super-predator population remains constant.

In [19]:
initial_populations_super = [
    100,    # prey
    50,     # predators
    10,     # super-predators
]

params_super = (
    1.0,    # a = prey_reproduction
    0.02,   # b = prey_predator_victim
    1.0,    # c = predator_mortality
    0.01,   # d = predator_reproduction
    0.01,   # e = prey_spredator_victim
    0.01,   # f = predator_spredator_victim
    0.0,    # g = spredator_mortality
    0.0,    # h = spredator_prey_reproduction
    0.0     # i = spredator_predator_reproduction
)

detail_plot_super(initial_populations=initial_populations_super, params=params_super)

Unstable Oscillation - Super Predators

By varying the parameters associated with the super predators, I can see how much their population affects the instability of the other populations in the system. This configuration I found particularly interesting, with plots showing how the super predator population is increasing with time, in turn decreasing the overall stability of the system.

In [20]:
times_super = np.linspace(0, 500, 5000)

initial_populations_super = [
    100,    # prey
    50,     # predators
    5,      # super-predators
]

params_super = (
    1.0,    # a = prey_reproduction
    0.02,   # b = prey_predator_victim
    1.0,    # c = predator_mortality
    0.01,   # d = predator_reproduction
    0.01,   # e = prey_spredator_victim
    0.01,   # f = predator_spredator_victim
    0.05,   # g = spredator_mortality
    0.0001, # h = spredator_prey_reproduction
    0.001   # i = spredator_predator_reproduction
)

detail_plot_super(initial_populations=initial_populations_super, params=params_super, times=times_super)

Unstable Oscillation - Predators Die Out

It was also possible to see that by decreasing the amount of prey deaths caused by super predators, this actually results in the exponential increase in super predator population beyond what the population of predators can sustain, resulting in the extinction of normal predators.

In [21]:
times_super = np.linspace(0, 600, 6000)

initial_populations_super = [
    100,    # prey
    50,     # predators
    5,      # super-predators
]

params_super = (
    1.0,    # a = prey_reproduction
    0.02,   # b = prey_predator_victim
    1.0,    # c = predator_mortality
    0.01,   # d = predator_reproduction
    0.001,  # e = prey_spredator_victim
    0.01,   # f = predator_spredator_victim
    0.05,   # g = spredator_mortality
    0.0001, # h = spredator_prey_reproduction
    0.001   # i = spredator_predator_reproduction
)

detail_plot_super(initial_populations=initial_populations_super, params=params_super, times=times_super)

Static Equilibrium

It is possible to solve the differential equations for the initial population, using those 3 parameters as independent variables. An example solution is shown below, where the differential equations all result in no change, leaving the system in static equilibrium.

In [24]:
times_super = np.linspace(0, 600, 6000)

initial_populations_super = [
    1000,    # prey
    100,     # predators
    10,      # super-predators
]

params_super = (
    1.100,   # a = prey_reproduction
    0.010,   # b = prey_predator_victim
    0.900,   # c = predator_mortality
    0.001,   # d = predator_reproduction
    0.010,   # e = prey_spredator_victim
    0.010,   # f = predator_spredator_victim
    1.100,   # g = spredator_mortality
    0.001,   # h = spredator_prey_reproduction
    0.001    # i = spredator_predator_reproduction
)

detail_plot_super(initial_populations=initial_populations_super, params=params_super, times=times_super)