# Python in Ecology – Stability Analysis of a Predator-Prey Model

As you may notice, my posts bounce back and forth between pop science, commentary, and programming/statistics/math. I basically blog on whatever I happen to be doing at the moment, or whatever I happen to be thinking about (but not everything I happen to be thinking about, which today involved alot of The NeverEnding Story because my dog kind of looks like Falcor, except cuter, and Spaceballs, because Mel Brooks).

Other that that, I’m working my way through this excellent book on modeling in biology by Sarah P. Otto and Troy Day: A Biologist’s Guide to Mathematical Modeling in Ecology and Evolution. Because I’m cheap/broke, I’m learning to use Python to work through some of the problems in this book, and I figured it might be good to share my experience for other interested ecologists who can’t afford Matlab/Mathematica/Maple or who embrace the open-source philosophy.

Today, I was working through one of their examples on equilibria and stability of predator-prey models (using a Type II functional response, albeit parameterized like a Monod function rather than the traditional Type II response equation). I thought I’d share how I got through the example using Python.

Here are the coupled, nonlinear ODE’s:

The first step is to define the variables and enter the equations. We define the variables as non-negative, although they can take on 0 values. Technically, r can be negative, but this won’t impact the results.

```import sympy as sm
R, P, r, K, c, e, a, d = sm.symbols('R, P, r, K, c, e, a, d', negative=False)

res = r*R*(1 - R/K) - P*c*R/(a+R)
cons = e*c*R*P/(a+R) - d*P

print res
print cons
```

The next step is to define the equalities, setting these equations equal to zero (which is the condition for equilibria).

```resEqual = sm.Eq(res, 0)
consEqual = sm.Eq(cons, 0)
```

Solving these ODEs individually would give us the equilibrium conditions for resources and consumers individually, but not simultaneously. The individual equilibria are null clines, where, for example, resources exhibit no change but predator abundances can still fluctuate. The simultaneous equilibria are the true equilibria, where neither resources nor consumers change. Finding the true equilibria is as simple as solving the system of equations simultaneously for both R and C.

```equilibria = sm.solve( (resEqual, consEqual), R, P )
print equilibria
```

This will return three equilibria, with two elements each. The elements are in the order we specified in sm.solve( ), R then P. So the first equilibrium, (0,0), is R = 0 and P = 0. This is intuitive. The system is at equilibrium when nothing is there. The third equilibrium, skipping the second for now, is (K/2 – a/2 + sqrt((K + a)^2)/2, 0). Sympy is a little too strict and absolute sometimes, but you can verify that this simplifies to (K,0), which is R = K and P = 0. When predators are absent, resources go to carrying capacity. The third equilibrium is (ad/(ce – d), -aer(-Kce + Kd + ad)/(K(ce – d)^2) ). This is the only equilibrium with both predators and resources present.

Next, we want to know some properties about these equilibria. Are they stable? Are they unstable? Since this model is purely symbolic, these questions are actually the same: What are the conditions of stability for these equilibria? We’ll use the Jacobian matrix to answer this question (I’m not covering what the Jacobian matrix is here, the book will tell you that).

First, set up a matrix containing the differential equations. Then, set up a matrix containing the variables you want to differentiate (in this case R and P). It is incredibly important that the equations and variables are entered in the same order, that is the resource equation goes first and the resource variable goes first.

```eqMat = sm.Matrix([ res, cons ])
resMat = sm.Matrix([ R, P ])

jacMat = eqMat.jacobian(resMat)
```

Next, we can substitute our equilibrium values into the Jacobian matrix and get the eigenvalues. These eigenvalues will tell us whether the conditions of stability. Start with the equilibrium R=0, P=0

```eq1mat = jacMat.subs([ (R,0), (P,0) ])
print eq1mat
```

This is a diagonal matrix with the elements r and -d, which means the eigenvalues are r and -d. d is the mortality rate of predators and cannot be negative. r is the population growth rate of resources and can be negative (if the population declines). Equilibrium stability requires that both eigenvalues be negative, meaning that populations will shrink towards the equilibrium. The only time the R=0,P=0 equilibrium is stable is when r < 0, or resources have a negative growth rate. Otherwise, resources will expand away from the origin at rate r.

Next, we examine the R=K, P=0 equilibrium.

```eq3mat = jacMat.subs([ (R,K), (P,0) ])
print eq3mat
```

This is an upper triangle matrix, meaning that the eigenvalues are again the diagonal elements, which are -r and Kce/(K + a) – d. Both eigenvalues are required to be negative for stability. So r must be > 1 (positive population growth rate), otherwise resources will fall away from the equilibrium if perturbed. Also, Kce/(K+a) – d must be < 0, which means Kce/(K+a) < d. The easiest way to interpret this is if c (the encounter rate between predators and prey) and e (the conversion efficiency of prey to predators once consumed) are sufficiently small. That is, if predators encounter prey rarely, and are very bad at making new predators, then predator populations will fail to establish. Otherwise, the equilibrium is unstable.

The final equilibrium is (ad/(ce – d), -aer(-Kce + Kd + ad)/(K(ce – d)^2) ), and this is the most complicated. Rather than type it all out, we’ll just substitute the values in from the saved equilibria.

```eq2mat = jacMat.subs([ (R, equilibria[1][0]), (P, equilibria[1][1]) ])

print eq2mat
```

This is a very complex matrix. Taking the eigenvalues results in a complex mess:

```evs = eq2mat.eigenvals()
print evs.keys()[0]
print evs.keys()[1]
```

NOTE: sympys saves the eigenvalues as a dictionary, so we have to print the dictionary keys to get the eigenvalues. That’s what the evs.keys( ) command is about.

These eigenvalues are incredibly complex, but they boil down to this:

-sqrt(some complicated function) / (2Kce(ce – d))

and

+sqrt(some complicated function) / (2Kce(ce – d))

We know two things. The square root value HAS to be positive. It CAN contain a -1, but that would just make the eigenvalue imaginary. If we assume the square root is positive, we can say that the eigenvalues are

-x / (2Kce(ce – d)) and +x / (2Kce(ce – d))

If the square root IS negative, then the eigenvalues are

-x * i /(2Kce(ce – d)) and +x * i / (2Kce(ce – d))

This makes no difference for stability (but it does make a difference whether the function oscillates or not). Regardless, notice that these eigenvalues have the same denominator, and are opposite in sign. The only way that + x / (2Kce(ce – d)) can be negative is if 2Kce(ce-d) is negative. However, if that were true, then -x / (2Kce(ce – d)) would be positive. Essentially, there’s no way both of these values can be negative, so this equilibrium is inherently unstable for any parameter values.

Let me know if I screwed that up, but I think it’s correct.