Logo

Analyzing more-dimensional systems

Example 1 - A Predator-Prey system

The following provides Python-code for analyzing the predator/prey system as discussed in Desch, Gertrud (2014/15): Einführung in die gleichungsbasierte Modellierung, chapter 7.2. (p. 41)

An ecosystem is the habitate of two species with one of them feeding on the other. Let $R(t)$ be the stock of predators and $B(t)$ be the stock of prey at time $t$. We assume:

  • both species have constant fertility $\beta_R, \beta_B$
  • natural mortality of both species (i.e. if it is independent of the other species) increases with high population densities: $\mu_R+\nu_R R(t), \mu_B+\nu_B B(t)$
  • the activity of predators increases the mortality of prey by $\kappa_B R(t)$
  • prey (as food source) reduces the mortality of predators by $\kappa_R B(t)$
  • without predators, the number of prey grows even in very small populations: $\delta_B:=\beta_B-\mu_B > 0$
  • without prey, the number of predators decreases: $\delta_R:=\mu_R-\beta_R > 0$

In terms of differential equations:

$$\frac{d}{dt}B(t) = (\beta_B-\mu_B) B(t) - \nu_B B(t)^2 - \kappa_B R(t) B(t)$$

$$:= B(t) (\delta_B - \nu_B B(t)- \kappa_B R(t))$$

$$\frac{d}{dt}R(t) = -(\mu_R-\beta_R) R(t) - \nu_R R(t)^2 + \kappa_R R(t) B(t))$$

$$:= R(t) (-\delta_R-\nu_R R(t)+\kappa_R B(t))$$

Parameters suggested are: $\delta_B=1, \delta_R=0.5,\kappa_B=0.1, \kappa_R=0.0025, \nu_B=0.0013, \nu_R=0.025$

The deployed Python-modules Matplotlib, Numpy, Scipy and Sympy are included in Python's Anaconda installation

1. Plot the system in time and in phase plane

iterating through different initial values

In [43]:
# import necessary libraries
import numpy as np
import sympy as sm
from scipy import integrate
import matplotlib.pyplot as plt
# show plots in notebook
% matplotlib inline

# prepare plots
fig = plt.figure(figsize=(15,5))
fig.subplots_adjust(wspace = 0.5, hspace = 0.3)
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)

# parameters
db = 1
dr = 0.5
kb = 0.1
kr = 0.0025
nb = 0.0013
nr = 0.025

def Sys(X, t=0):
    # here X[0] = B and X[1] = R    
    return np.array([ X[0] * (db - nb*X[0] - kb*X[1]) , X[1] * (-dr - nr*X[1] + kr*X[0]) ])

# generate 1000 linearly spaced numbers for x-axes
t = np.linspace(0, 50,  1000)

# initial values:
iv = np.linspace(0, 50, 10)

# plot for different initial values
for n,i in enumerate(iv):
    # define intial values
    Sys0 = np.array([iv[n], iv[-n]])
    
    # integrate
    # type "help(integrate.odeint)" if you want more information about integrate.odeint inputs and outputs.
    X, infodict = integrate.odeint(Sys, Sys0, t, full_output=True)
    # infodict['message']                      # integration successful
    # transpose output
    x,y = X.T
    # plot
    ax1.plot(x, 'b-')
    ax1.plot(y, 'r-')
    ax2.plot(x, y, 'k')

ax1.plot(0, 'b-', label='Prey')
ax1.plot(0, 'r-', label='Predator')
ax1.set_title("Dynamics in time")
ax1.set_xlabel("time")
ax1.grid()
ax1.legend(loc='best')

ax2.set_xlabel("Prey")
ax2.set_ylabel('Predator')  
ax2.set_title("Phase space")
ax2.grid()

# aproximate equilibrium
print(x[-1], y[-1])
(265.48550383988601, 6.5486583553093283)

2. Find eqilibria (aka fixed points aka extrema)

In [44]:
# define the system symbolically (assuming a predator-prey-system with no negative values)
r, b = sm.symbols('r, b', negative=False)

# differential equations
B = b * (db - nb*b - kb*r)
R = r * (-dr - nr*r + kr*b)

# use sympy's way of setting equations to zero
BEqual = sm.Eq(B, 0)
REqual = sm.Eq(R, 0)

# compute equilibria (fixed points, extrema)
fp = sm.solve( (BEqual, REqual), b, r )
print('The equilibria of this system are %s' %fp)
The equilibria of this system are [(0.0, 0.0), (265.486725663717, 6.54867256637168), (769.230769230769, 0.0)]

3. Draw nullclines

In [45]:
# solve zero-set-equations independently from each other
bcline = sm.solve(BEqual, b)
rcline = sm.solve(REqual, r )
# take only non-zero solution (needs to be casted into string)
bc = str([u for u in bcline if u != 0.])
rc = str([v for v in rcline if v != 0.])

# solve bc for y (= the predator values)
Bcline = []
for r in y:
    # needs to be evaluated from string
    Bcline.append(eval(bc))

# solve rc for x (= the prey values)
Rcline = []
for b in x:
    # needs to be evaluated from string
    Rcline.append(eval(rc))

#prepare plot
fig = plt.figure(figsize=(8,6))
ax3 = fig.add_subplot(1,1,1)
# plot B-nullcline on x-axes   
ax3.plot(Bcline, y, 'r-', label = 'Prey-nullcline')
# plot R-nullcline on y-axes   
ax3.plot(x, Rcline, 'b-', label = 'Predator-nullcline')
ax3.set_title("Phase space with nullclines")
ax3.set_xlabel("Prey")
ax3.set_ylabel("Predator")
ax3.plot(x, y, 'k-')
ax3.grid()
# plot fixed points 
for f in fp:
    ax3.plot(f[0], f[1],"red", marker = "o", markersize = 10.0)
ax3.plot(0, 0, 'ro', label = 'Equilibria')
ax3.legend(loc='best')
Out[45]:
<matplotlib.legend.Legend at 0x11598dd0>

4. Draw a stream-plot (aka quiver-plot)

In [46]:
fig2 = plt.figure(figsize=(8,6))
ax4 = fig2.add_subplot(1,1,1)

# define dimensions for grid
u = np.linspace(0, 800, 30)
v = np.linspace(-15, 25, 20)
# create grid
X1 , Y1  = np.meshgrid(u, v)
# evaluate system on the grid
DX1, DY1 = Sys([X1, Y1])
# norm results
M = (np.hypot(DX1, DY1))
# avoid zero division errors
M[M == 0] = 1.
# normalize each arrow
DX1 /= M
DY1 /= M

# plot stream-plot
ax4.quiver(X1, Y1, DX1, DY1, M, pivot='mid')

# same as above
ax4.plot(Bcline, y, 'r-', label = 'Prey-nullcline')
ax4.plot(x, Rcline, 'b-', label = 'Predator-nullcline')
ax4.set_title("Quiverplot with nullclines")
ax4.set_xlabel("Prey")
ax4.set_ylabel("Predator")
ax4.plot(x, y, 'k-')
ax4.grid()
ax4.legend(loc='best')
for f in fp:
    ax4.plot(f[0], f[1],"red", marker = "o", markersize = 10.0)
ax4.plot(0, 0, 'ro', label = 'Equilibria')
ax4.legend(loc='best')
Out[46]:
<matplotlib.legend.Legend at 0x1179f4b0>

5. Calculate eigenvalues from Jacobian matrix

$A = \begin{pmatrix} \frac{df}{dx} & \frac{df}{dy} \\ \frac{dg}{dx} & \frac{dg}{dy} \\ \end{pmatrix} = \begin{pmatrix} a_{11} & a_{12} \\ a_{21} & a_{22} \\ \end{pmatrix}$

Determinant: $detA = a_{11}a_{22} - a_{12}a_{21}$

Trace: $trA=a_{11} + a_{22}$

Eigenvalues: $\lambda = \frac{trA \pm \sqrt{(trA)^2-4detA}}{2}$

In [47]:
# repeat symbol definition (see step 2 above)
r, b = sm.symbols('r, b', negative=False)

# put equations into matrix-form and compute Jacobian-matrix  
eqMat = sm.Matrix([ B, R ])
Mat = sm.Matrix([ b, r ])
jacMat = eqMat.jacobian(Mat)
print('Jacobian matrix: \n %s' % np.matrix(jacMat))
Jacobian matrix: 
 [[-0.0026*b - 0.1*r + 1 -0.1*b]
 [0.0025*r 0.0025*b - 0.05*r - 0.5]]
In [48]:
# needed for eigenvalue computation
import numpy.linalg as linalg

# iterate through equilibria (calulated in step 2 above)
for item in fp:
    # substitute equilibrium values into the Jacobian matrix and get the eigenvalues
    # needs to be converted to float, with .astype(np.float64)
    eqmat = np.array(jacMat.subs([ (b, item[0]), (r, item[1])])).astype(np.float64)
    # eigenvalues as first part of linalg.eig()
    eigenValues,eigenVectors = linalg.eig(eqmat)
    # take only real parts, with eigenValues.real
    print('The real parts of the eigenvalues for the fixed point (%s, %s) are: %s' 
          %(item[0], item[1], eigenValues.real))
    if eigenValues.real[0] < 0 and eigenValues.real[1] < 0:
        print('The fixed point is a sink, the equilibrium is stable')
    if eigenValues.real[0] > 0 or eigenValues.real[1] > 0:
        print('The fixed point is a source, the equilibrium is unstable')
    if eigenValues.real[0] > 0 and eigenValues.real[1] < 0:
        print('The fixed point is a saddle')
    # find dominant one (the largest absolut)
    domEV = max([abs(x) for x in eigenValues.real])
    # find its position
    posdomEV = [abs(x) for x in eigenValues.real].index(domEV)
    print('The dominant eigenvalue is the %s one:' %(posdomEV + 1))
    print('-------------------------------------------')
The real parts of the eigenvalues for the fixed point (0.0, 0.0) are: [ 1.  -0.5]
The fixed point is a source, the equilibrium is unstable
The fixed point is a saddle
The dominant eigenvalue is the 1 one:
-------------------------------------------
The real parts of the eigenvalues for the fixed point (265.486725663717, 6.54867256637168) are: [-0.25442478 -0.25442478]
The fixed point is a sink, the equilibrium is stable
The dominant eigenvalue is the 1 one:
-------------------------------------------
The real parts of the eigenvalues for the fixed point (769.230769230769, 0.0) are: [-1.          1.42307692]
The fixed point is a source, the equilibrium is unstable
The dominant eigenvalue is the 2 one:
-------------------------------------------

Example 2 - An enzymatic reaction

The following provides Python-code for analyzing a model for the enzymatic reaction $S+E <-> K -> S+P$ as discussed in Desch, Gertrud (2014/15): Einführung in die gleichungsbasierte Modellierung chapter 7.2.7 p.46

The model assumes an enzyme $E$ in solution. A substrate $S$ is constantly generated and enzymatically degraded according to the following reaction:

$S+E <-> K -> S+P$

$S$ and $E$ form a complex $K$, which in turn can decompose either into $S$ and $E$, or into $P+Q+E$, breaking the substrate down into the two products $P$ and $Q$. Products $P$ and $Q$ precipitate from solution and do not take part in any further reactions. The enzyme bound in the complex is now free for new reactions.

The concentration of the substrate (in mol/mL) is referred to by $c_S$, $c_E$ is the concentration of the free enzyme, and $c_K$ is the concentration of the complex. $c_G$ is the concentration of all free and bound enzymes, i.e. $c_G = c_E + c_K$. This quantity is constant.

Every second and milliliter solution, $b_S mol$ of the substrate are freshly generated. According to the law of mass action, $k_1c_Sc_E mol$ of the complex are formed every second and milliliter of solution, and the corresponding amount of substrate and free enzyme are bound. Per second and milliliter of solution, $k_2c_K mol$ of the complex decompose to $S+E$, and $k_3c_K mol$ of the complex decompose to $S+P+Q$. We set up a model for the enzyme and substrate concentrations.

We refer to one milliliter of solution. Flows are:

  • New formation of substrate: $b_S \frac{mol}{mL s}$ – is added to $c_S$.
  • Complex formation: $k_1 c_Sc_E \frac{mol}{mL s}$ – is added to $c_K$ and substracted from $c_E$ and $c_S$.
  • Complex decomposition into $S+E: k_2 c_K \frac{mol}{mL s}$ – is removed from $c_K$ and added to $c_E$ and $c_S$.
  • Complex decomposition into $P+Q+E: k_3 c_K \frac{mol}{mL s}$ – is removed from $c_K$ and added to $c_E$.

$\frac{d}{dt}c_S(t) = b_s - k_1 c_S(t) c_E(t)+k_2 c_K(t)$

$\frac{d}{dt}c_E(t) = -k_1 c_s(t) c_E(t)+(k_2+k_3) c_K(t)$

$\frac{d}{dt}c_K(t) = k_1 c_S (t)c_E (t)- (k_2 + k_3)c_K(t)$

From this, the differential equaion for $c_G = c_E + c_K$ is obtained

$\frac{d}{dt} c_G(t) = \frac{d}{dt} c_E(t) + \frac{d}{dt} c_K (t) = 0$

Hence, $c_G$ is constant, and the term $c_K(t)$ in the system can be replaced by $c_G-c_E(t)$. As a result, a system with two unknown functions remains:

$\frac{d}{dt}c_S(t) = b_s - k_1*c_S(t)*c_E(t)+k_2*(c_G-c_E(t))$

$\frac{d}{dt}c_E(t) = -k_1*c_s(t)*c_E(t)+(k_2+k_3)*(c_G-c_E)(t))$

1. Plot the system in time and in phase plane

In [33]:
# import necessary libraries
import matplotlib.pyplot as plt
# show plots in notebook
% matplotlib inline
import numpy as np
import sympy as sm
from scipy import integrate

# prepare plots
fig = plt.figure(figsize=(15,5))
fig.subplots_adjust(wspace = 0.5, hspace = 0.3)
ax1 = fig.add_subplot(1,2,1)
ax2 = fig.add_subplot(1,2,2)

bs = 0.1
k1 = 0.4
k2 = 0.5
k3 = 0.2
cg = 1
# initial values:
S = 0.5
E = 0.9

def Sys(X, t=0):
    # here X[0] = x and x[1] = y    
    return np.array([ bs - k1 * X[0] * X[1] + k2 * (cg - X[1]) , - k1 * X[0] * X[1] + (k2 + k3) * (cg - X[1]) ])

# generate 1000 linearly spaced numbers for x-axes
t = np.linspace(0, 200,  100)

Sys0 = [S, E]

# type "help(integrate.odeint)" if you want more information about integrate.odeint inputs and outputs.
X, infodict = integrate.odeint(Sys, Sys0, t, full_output=True)
# infodict['message']                      # integration successful

x,y = X.T

ax1.plot(x, 'b-', label='substrat')
ax1.plot(y, 'r-', label='enzyme')
ax2.plot(x, y, 'ko')

ax1.set_title("Dynamics in time")
ax1.set_xlabel("time")
ax1.grid()
ax1.legend(loc='best')

ax2.set_xlabel("substrat")
ax2.set_ylabel('enzyme')  
ax2.set_title("Phase space portrait")
ax2.grid()

# approximate equilibria
print x[-1]
print y[-1]
1.74383477076
0.500898688747

2. Find eqilibria (aka fixed points)

In [37]:
# equilibria with sympy
import sympy as sm

# define parameter-symbols
s, e = sm.symbols('s, e', negative=False)

# equations
S = bs - k1 * s * e + k2 * (cg - e)
E = - k1 * s * e + (k2 + k3) * (cg - e)

# use sympy's way of setting equations to zero
SEqual = sm.Eq(S, 0)
EEqual = sm.Eq(E, 0)

# compute fixed points
equilibria = sm.solve( (SEqual, EEqual), s, e )
print('The fixed point(s) of this system are: %s' %equilibria)
The fixed point(s) of this system are: [(1.75000000000000, 0.500000000000000)]

3. Draw nullclines and quiver-plot

In [35]:
# find nullclines 

# parameter-symbols have to be defined again (if re-compiling of code above is to be evaded)
s, e = sm.symbols('s, e', negative=False)

# solve zero-set-equations independently from each other
scline = sm.solve(SEqual, s)
sc = str(scline)

ecline = sm.solve(EEqual, e )
ec = str(ecline)

# solve bc for y (= the predator values)
Scline = []
for e in y:
    # needs to be evaluated from string, comes as a list -> take only first item
    Scline.append(eval(sc)[0])

# solve rc for x (= the prey values)
Ecline = []
for s in x:
    Ecline.append(eval(ec)[0])

fig = plt.figure(figsize=(8,6))
ax3 = fig.add_subplot(1,1,1)
# plot B-nullcline on x-axes   
ax3.plot(Scline, y, 'b-', label = 'Substrat-nullcline')
# plot R-nullcline on y-axes   
ax3.plot(x, Ecline, 'r-', label = 'Enzyme-nullcline')
ax3.set_title("Phase space with nullclines and Quiverplot")
ax3.set_xlabel("Substrat")
ax3.set_ylabel("Enzyme")
ax3.plot(x, y, 'ko', label = 'Phase space portrait')
ax3.grid()
for f in equilibria:
    fpx = f[0] 
    fpy = f[1]
    ax3.plot(fpx, fpy,"red", marker = "o", markersize = 10.0)
ax3.plot(fpx, fpy, 'ro', label = 'Equilibria')
ax3.legend(loc='best')

# quiverplot
# define a grid and compute direction at each point
# use start and end of all possible values
u = np.linspace(min(float(fpx), Ecline[0], Scline[0], x[0]), max(float(fpx), Scline[-1], x[0], x[-1]), 30)
v = np.linspace(min(float(fpy), Ecline[0], Ecline[-1], y[0]), max(float(fpy), Ecline[0], Scline[0], y[0]), 20)

#print(float(fpx),float(fpy), Scline[0], Ecline[0], Scline[-1], Ecline[-1], y[0], x[0], y[-1], x[-1])

X1 , Y1  = np.meshgrid(u, v)                    # create a grid
DX1, DY1 = Sys([X1, Y1])                        # compute growth rate on the grid
M = (np.hypot(DX1, DY1))                        # norm growth rate 
M[ M == 0] = 1.                                 # avoid zero division errors 
DX1 /= M                                        # normalize each arrows
DY1 /= M

ax3.quiver(X1, Y1, DX1, DY1, M, pivot='mid')
Out[35]:
<matplotlib.quiver.Quiver at 0x118e3dd0>

4. Check stability of Equilibria

In [38]:
# put equations into matrix-form and compute Jacobian-matrix  
eqMat = sm.Matrix([ S, E ])
Mat = sm.Matrix([ s, e ])
jacMat = eqMat.jacobian(Mat)
print('Jacobian Matrix: \n %s' % np.matrix(jacMat))

import numpy.linalg as linalg

for item in equilibria:
    # substitute equilibrium values into the Jacobian matrix and get the eigenvalues
    # needs to be converted to float, with .astype(np.float64)
    eqmat = np.array(jacMat.subs([ (s, item[0]), (e, item[1])])).astype(np.float64)
    # eigenvalues as first part of linalg.eig()
    eigenValues,eigenVectors = linalg.eig(eqmat)
    # only the real parts
    print('The real parts of the eigenvalues for the fixed point (%s, %s) are %s:' 
          %(item[0], item[1], eigenValues.real))
    if eigenValues.real[0] < 0 and eigenValues.real[1] < 0:
        print('The fixed point is a sink, the equilibrium is stable')
    if eigenValues.real[0] > 0 or eigenValues.real[1] > 0:
        print('The fixed point is a source, the equilibrium is unstable')
    if eigenValues.real[0] > 0 and eigenValues.real[1] < 0:
        print('The fixed point is a saddle')
    # find dominant one (the largest absolut)
    domEV = max([abs(x) for x in eigenValues.real])
    posdomEV = [abs(x) for x in eigenValues.real].index(domEV)
    print('The dominant eigenvalue is the %s one:' %(posdomEV + 1))
    print('-------------------------------------------')
Jacobian Matrix: 
 [[-0.4*e -0.4*s - 0.5]
 [-0.4*e -0.4*s - 0.7]]
The real parts of the eigenvalues for the fixed point (1.75000000000000, 0.500000000000000) are [-0.02540333 -1.57459667]:
The fixed point is a sink, the equilibrium is stable
The dominant eigenvalue is the 2 one:
-------------------------------------------

5. Iterate through a range of $k_3$ and see what the dominant Eigenvalue does

In [42]:
# prepare plots
fig = plt.figure(figsize=(15,5))
fig.subplots_adjust(wspace = 0.8, hspace = 0.3)
ax4 = fig.add_subplot(1,2,1)
ax5 = fig.add_subplot(1,2,2)

K = np.linspace(0.2, 0.9, 10)

dEV = []

for k3 in K:
    X, infodict = integrate.odeint(Sys, Sys0, t, full_output=True)
    x,y = X.T
    ax4.plot(x, y, label = 'k3 = %s' %k3)
    
    # define parameter-symbols for solution
    s, e = sm.symbols('s, e', negative=False)

    # equations
    S = bs - k1 * s * e + k2 * (cg - e)
    E = - k1 * s * e + (k2 + k3) * (cg - e)

    # use sympy's way of setting equations to zero
    SEqual = sm.Eq(S, 0)
    EEqual = sm.Eq(E, 0)

    # compute fixed points
    equilibria = sm.solve( (SEqual, EEqual), s, e )
    #print equilibria

    eqMat = sm.Matrix([ S, E ])
    Mat = sm.Matrix([ s, e ])
    jacMat = eqMat.jacobian(Mat)
    #print('Jacobian Matrix: \n %s' % np.matrix(jacMat))

    # substitute equilibrium values into the Jacobian matrix and get the eigenvalues
    # make sure equilibria exist
    try:
        # needs to be converted to float, with .astype(np.float64)
        eqmat = np.array(jacMat.subs([ (s, equilibria[0][0]), (e, equilibria[0][1])])).astype(np.float64)
        # eigenvalues as first part of linalg.eig()
        eigenValues,eigenVectors = linalg.eig(eqmat)
        # find dominant one (the largest absolut)
        domEV = max([abs(x) for x in eigenValues.real])
        # get its position
        posdomEV = [abs(x) for x in eigenValues.real].index(domEV)
        ax5.plot(k3, eigenValues.real[posdomEV],"red", marker = "o", markersize = 10.0)
    except:
        pass

ax4.set_title('Development of phase space portrait')
ax4.legend(loc='center left', bbox_to_anchor=(1, 0.5))
ax5.set_title('Development of dominant eigenvalues')
Out[42]:
<matplotlib.text.Text at 0x17210a30>

Example 3 - A 3-dimensional predator-prey-system

The following provides Python-code for analyzing the model as discussed in: Van Kooten et al. (2005). Bistability and an Allee effect as emergent consequences of stage-specific predation. Journal of Theoretical Biology 237, 67–74.

This model analyzes a so called emergent Allee-effect, an Allee-effect that is not an a priori assumption of the model (like here), but results from an interaction between a size- or stage-dependent predator and density-dependent development of its prey population.

In order to account for this interaction, the model distinguishes between juvenile and adult prey individuals and considers predators that feed on adult prey. This is expressed in the following three ordinary differential equations for the development of, respectively, juvenile (J), adult (A), and predator density (P):

$\frac{dJ}{dt} = \beta A - \frac{J}{1+J^2}-\mu_J J$

$\frac{dA}{dt} = \frac{J}{1+J^2} -\mu_A A - A P$

$\frac{dP}{dt} = P(\epsilon A - \mu_P)$

The term $\frac{J}{1+J^2}$ indicates the maturation rate of the juvenile population, implying that maturation goes to zero as juvenile density goes to infinity. The underlying energy budget model assumes that a fixed fraction of ingested energy is spent on maintenance plus growth. In cases of low food densities ingested energy just suffices to cover maintenance requirements. Maturation slows down or might even stop.

Originally, the model refers to the observation that predators who feed specifically on the small individuals of a size-structured population act in a way like gardeners weeding out most of the small individuals so that the remaining ones can grow faster and produce more offspring which serves as prey. It has been shown however, that the Allee-effect can occur no matter which fraction of the size-structured population is reduced. In this case therefore, the model considers predators that feed on adult prey.

Suggested parameters are:

  • Generation rate of Juveniles per capita of adult: $\beta = 1.0$
  • Reproduction rate of predators per capita of adult: $\epsilon=1.0$
  • Mortality rate of juveniles: $\mu_J = 0.05$
  • Mortality rate of adults: $\mu_A = 0.1$
  • Mortality rate of predators: $\mu_P = 0.56$
In [1]:
# import libraries
import numpy as np
from scipy import integrate
import sympy as sm
import matplotlib.pyplot as plt
# show plots in notebook
% matplotlib inline

beta = 1.
eps = 1.
muj = 0.05
mua = 0.1

def Sys(X, t=0):
    # here X[0] = J, x[1] = A, and X[2] = P     
    return np.array([ beta*X[1] - X[0]/(1 + X[0]**2) - muj*X[0] , 
                     X[0]/(1+X[0]**2) - X[1]*X[2] - mua*X[1], 
                     X[1]*X[2]*eps - mup*X[2] ])

# generate 1000 linearly spaced numbers for x-axes
t = np.linspace(0, 50,  200)

# initial values
Sys0 = [1, 1, 1] # = [J, A, P]

# prepare plots
fig = plt.figure(figsize=(15,5))
fig.subplots_adjust(wspace = 0.5, hspace = 0.3)
ax = fig.add_subplot(1,2,1)
ax0 = fig.add_subplot(1,2,2)

# generate linearly spaced numbers for iterating mup
MM = np.linspace(0.4, 0.6, 2)

for mup in MM:
    X, infodict = integrate.odeint(Sys, Sys0, t, full_output=True)    
    j,a,p = X.T

    ax.plot(j, 'g-')
    ax.plot(a, 'b-')
    ax.plot(p, 'r-')

    ax0.plot(j, a, '.', label='$\mu_p = $ %s' %mup)
    
ax.plot(0, 'g-', label='Juvenile')
ax.plot(0, 'b-', label='Adult')
ax.plot(0, 'r-', label='Predator')

ax.set_title('Dynamics in time for $\mu_p =$ %s' %MM)
ax.set_xlabel('time')
ax.grid()
ax.legend(loc='best')
ax0.legend(loc='best')

ax0.set_xlabel('Juvenile')
ax0.set_ylabel('Adult')  
ax0.set_title('Phase space portrait')
ax0.grid()

#approximate equilibria
#print j[-1]
#print a[-1]
#print p[-1]

Note that the population of juveniles seems to develope very differently, depending on whether $\mu_p = 0.4$ or $\mu_p = 0.6$

Lets check the equilibrium for $\mu_p = 0.5$

2. Find equilibria

In [12]:
beta = 1
eps = 1
muj = 0.05
mua = 0.1
mup = 0.5

#define the variables as non-negative and exclude non-real results
j, a, p  = sm.symbols('j, a, p', negative=False, real = True)
J = beta * a - j / (1 + j**2) - muj * j
A = j / (1 + j**2) - a * p - mua * a
P = a * p * eps - mup * p
#print J
#print A
#print P

# use sympy's way of setting equations to zero
JEqual = sm.Eq(J, 0)
AEqual = sm.Eq(A, 0)
PEqual = sm.Eq(P, 0)

#equilibria
equilibria = sm.solve((JEqual, AEqual, PEqual), j, a, p)
print equilibria
[(0.0, 0.0, 0.0), (0.683375209644600, 0.500000000000000, 0.831662479035540), (2.00000000000000, 0.500000000000000, 0.700000000000000), (7.31662479035540, 0.500000000000000, 0.168337520964460), (13.3790881602597, 0.743282675569981, 0.0)]

Note that there are five equilibria with $\mu_p = 0.5$ So lets check for the stability of these euqilibria.

3. Check stability of equilibria

In [13]:
# put equations into matrix-form and compute Jacobian-matrix  
eqMat = sm.Matrix([ J, A, P ])
Mat = sm.Matrix([ j, a, p ])
jacMat = eqMat.jacobian(Mat)
print('Jacobian Matrix \n %s' % np.matrix(jacMat))

# needed for computing eigenvalues
import numpy.linalg as linalg

for item in equilibria:
    # substitute equilibrium values into the Jacobian matrix and get the eigenvalues
    # needs to be converted to float, with .astype(np.float64)
    eqmat = np.array(jacMat.subs([ (j, item[0]), (a, item[1]), (p, item[2])])).astype(np.float64)
        
    # eigenvalues as first part of linalg.eig()
    eigenValues,eigenVectors = linalg.eig(eqmat)
    print('-------------------------------------------')
    # take only real parts
    print('The real parts of the eigenvalues for the fixed point (%s, %s, %s) are %s:' 
          %(item[0], item[1], item[2], eigenValues.real))
    # check sign of real parts
    if all([d < 0 for d in eigenValues.real]): 
        print('This fixed point is a sink, the equilibrium is stable')
    else:
        print('The fixed point is a source, the equilibrium is unstable')
    # find dominant eigenvalue (the largest absolut)
    domEV = max([abs(x) for x in eigenValues.real])
    posdomEV = [abs(x) for x in eigenValues.real].index(domEV)
    print('The dominant eigenvalue is the %s one:' %(posdomEV + 1))
Jacobian Matrix 
 [[2*j**2/(j**2 + 1)**2 - 0.05 - 1/(j**2 + 1) 1 0]
 [-2*j**2/(j**2 + 1)**2 + 1/(j**2 + 1) -p - 0.1 -a]
 [0 p a - 0.5]]
-------------------------------------------
The real parts of the eigenvalues for the fixed point (0.0, 0.0, 0.0) are [-1.68207949  0.53207949 -0.5       ]:
The fixed point is a source, the equilibrium is unstable
The dominant eigenvalue is the 1 one:
-------------------------------------------
The real parts of the eigenvalues for the fixed point (0.683375209644600, 0.500000000000000, 0.831662479035540) are [-0.88371406 -0.1728067  -0.1728067 ]:
This fixed point is a sink, the equilibrium is stable
The dominant eigenvalue is the 1 one:
-------------------------------------------
The real parts of the eigenvalues for the fixed point (2.00000000000000, 0.500000000000000, 0.700000000000000) are [ 0.05371647 -0.39185823 -0.39185823]:
The fixed point is a source, the equilibrium is unstable
The dominant eigenvalue is the 2 one:
-------------------------------------------
The real parts of the eigenvalues for the fixed point (7.31662479035540, 0.500000000000000, 0.168337520964460) are [-0.02635119 -0.13716067 -0.13716067]:
This fixed point is a sink, the equilibrium is stable
The dominant eigenvalue is the 2 one:
-------------------------------------------
The real parts of the eigenvalues for the fixed point (13.3790881602597, 0.743282675569981, 0.0) are [-0.07225309 -0.07225309  0.24328268]:
The fixed point is a source, the equilibrium is unstable
The dominant eigenvalue is the 3 one:

The first and the last equilibria are less interesting. At least one of the populations is extinct in them. The second, third and forth equilibrium however seem interesting, with the middle one being unstable.

Lets check a range of $\mu_P$-values to find out more about these equilibria.

4. Iterate through different values for $\mu_P$ and check stability of equilibria

We consider only those with all populations alive.

In [4]:
Eq = []

#create an array of mup-values
M = np.arange(0.35, 0.65, 0.01)

for mup in M:
    #define the variables as non-negative
    j, a, p  = sm.symbols('j, a, p', negative=False)
    J = beta * a - j / (1 + j**2) - muj * j
    A = j / (1 + j**2) - a * p - mua * a
    P = a * p * eps - mup * p

    # use sympy's way of setting equations to zero
    JEqual = sm.Eq(J, 0)
    AEqual = sm.Eq(A, 0)
    PEqual = sm.Eq(P, 0)

    # get equilibria
    equilibria = sm.solve((JEqual, AEqual, PEqual), j, a, p)

    # get rid of those with one of the species being zero
    equili = []
    for ii in equilibria:
        if not any([i == 0.0 for i in ii]):
            equili.append(ii)

    # in some cases sm.solve produces a formate which has to be converted into the image of a complex number, 
    # consisting just of real-part and complex part 
    RE = []
    for item in equili:
        reqi = [f.as_real_imag() for f in item]
        re = []
        for it in reqi:
            # here we take out the real-part
            re.append(it[0])
        RE.append(re)        
    Eq.append(RE)    

Note that this produces some duplicated non-real equilibria with no significance. We have to get rid of these.

In [5]:
# get rid of duplicates in Eq (as a function)
def remove_duplicates(List):
    L = []
    for item in List: 
        LL = []
        for ii in item:
            # append only if its not there
            if ii not in LL:
                LL.append(ii)
            # if it should be there already, remove it
            else:
                LL.remove(ii)
        L.append(LL)
    return L

Equ = remove_duplicates(Eq)

Separate equilibria into stable and unstable parts

In [31]:
# http://matplotlib.org/users/gridspec.html
import matplotlib.gridspec as gridspec
# prepare plots
fig = plt.figure(figsize = (18, 12))
gs = gridspec.GridSpec(3, 3)
plot_coor = [[1,0],[1,1],[1,2]]
#ax3 = fig.add_subplot(2,2,2)
plot_coor2 = [[2,0],[2,1],[2,2]]

population = ['Juvenile','Adult','Predator']

for q, pop in enumerate(population):    
    stabEq = []
    unstabEq = []
    domEigV = []
    xs = []
    xu = []
    xd = []

    for n,item in enumerate(Equ):
        for i in item: 
            eqmat = np.array(jacMat.subs([ (j, i[0]), (a, i[1]), (p, i[2])])).astype(np.float64)
            # eigenvalues as first part of linalg.eig()
            eigenValues,eigenVectors = linalg.eig(eqmat)
            # check stability
            if all([d < 0 for d in eigenValues.real]): 
                # stable
                stabEq.append(i[q])
                xs.append(M[n])
            else:
                # unstable
                unstabEq.append(i[q])
                xu.append(M[n])
            # find dominant one (the largest absolut)
            domEV = max([abs(x) for x in eigenValues.real])
            # get its position
            # posdomEV = [abs(x) for x in eigenValues.real].index(domEV)
            domEigV.append(eigenValues.real[0]) #posdomEV])
            xd.append(M[n])

    ax2 = plt.subplot(gs[plot_coor[q][0],plot_coor[q][1]])
    ax2.plot(xs, stabEq, 'ro', label = 'stable')
    ax2.plot(xu, unstabEq, 'bo', label = 'unstable')
    ax2.grid()
    ax2.legend(loc='best')
    ax2.set_xlabel('$\mu_P$')
    ax2.set_ylabel(pop)  
    ax2.set_title('Equilibrium of %s in respect to $\mu_P$' %pop)
ax3 = plt.subplot(gs[plot_coor2[0][0],plot_coor2[0][1]])
ax3.plot(xd, domEigV,'ro')
ax3.grid()
ax3.set_xlabel('$\mu_P$')
ax3.set_title('Development of dominant Eigenvalue')
    
Out[31]:
<matplotlib.text.Text at 0x15ab0770>