Predator-Prey Models

Introduction

Volterra (1926) proposed a simple two-component model to explain oscillatory populations of certain Adriatic fish.  Similar equations were also proposed by Lotka (1920) as a model of chemical oscillations.  Let X represent the population of a prey fish and Y represent the population of a predator.  The Lotka-Volterra model takes the form

d(X/d)τ == (A - B Y) 
        X
d(Y/d)τ == (C X - D) 
        Y

where A represents the growth of the prey population in the absence of predators, B represents the population loss due to predation, D represents the death rate of predators in the absence of prey, C represents the growth of the predator population, and τ is a time variable.  The interaction between these species is taken to be proportional to both populations.  The model is obviously too simple — finite environmental resources would not permit exponential growth of the prey population even if predators were absent and there must be some saturation mechanism for the predator population — but it is a useful starting point for predator-prey models.

This is an example of an autonomous dynamical system that is described by a system of first-order ordinary differential equations of the form

(∂_t) r  =  v[r]

where r = {(Subscript[r, 1]), (Subscript[r, 2,]) ⋯, 
          (Subscript[r, n])} is point in an n-dimensional phase space and v = {(Subscript[v, 1]), (Subscript[v, 2]), ⋯, (Subscript[v, n])} is an n-dimensional velocity field that depends only upon the coordinates r and has no explicit dependence upon time.  Thus, if one knows the phase-space coordinates at any time the subsequent evolution of the system can, at least in principle, be computed for all later times.  The dynamics of the system are described by a family of nonintersecting trajectories, or orbits, in phase space.  The Lotka-Volterra equations are a two-dimensional autonomous system describing the phase-space trajectory r[τ] = {X[τ], Y[τ]}.

Points (Subscript[r, fp]) where the velocity vanishes, such that v[(Subscript[r, fp])] == 0, are described as fixed points and are especially important in the analysis of autonomous systems.  Near any fixed point, we can linearize the equations of motion by performing a series expansion of the velocity field.  The resulting linearized equations of motion can then be expressed in the form

r[t]  ≈  
      (Subscript[r, fp])  ⟹   (∂_t) r[
            t]  ≈  
      S·(r[t] - (Subscript[r, fp]))

where

\(TraditionalForm\`S\_\(i, j\) = \((∂v\_i\/∂r\_j)\)\_fp\)

is the stability matrix and is composed of derivatives of velocity with respect to position evaluated at the fixed point.  It is useful to define a displacement vector

\(TraditionalForm\`δ[t]\  = 
      r[t] - r\_fp\ \ \ ⟹\ \ \ \(∂\_t δ[t]\)\  == \ S\ δ[t]\)

and a unitary matrix U that diagonalizes S, such that

\(TraditionalForm\`U\ S\ U\^\(-1\)\  = \ λ\ , \ \ η\  = \ 
      U\ δ\ \ ⟹\ \ \ \ \(∂\_t η\_i[t]\)\  == \ λ\_i\ η\_i[t]\)

permits explicit solutions in the form

(Subscript[η, i])[
        t]  =  (Subscript[η, i])[0] Exp[(Subscript[λ, i]) t]

where (Subscript[λ, i]) is an eigenvalue of S and (Subscript[η, i]) is the corresponding eigenvector.  The dynamical behavior of the system near its fixed points can now be classified according to the properties of the eigenvalues of the stability matrix.  For simplicity, consider a two-dimensional system such that the eigenvalues are either both real or are a pair of complex conjugates.  Although one can classify all possibilities rather completely, for the present purposes it is sufficient to consider the case of complex eigenvalues.  If the real part of an eigenvalue is negative, the corresponding fixed point is considered stable because small initial displacements are damped exponentially, whereas if the real part is positive the fixed point is unstable because small initial displacements grow exponentially.  The imaginary parts of these eigenvalues produce rotations in phase space, either clockwise or counterclockwise spirals.  When Re[λ] = 0, the orbit reduces to a circle at fixed distance that is described as neutrally stable.

In this notebook we study some of the behavior of predator-prey models.  We begin with the original Lotka-Volterra equations and use dimensional analysis to reduce the number of parameters from 4 to 1 by rescaling the population and time variables.  We then determine the fixed points, analyze their stability, and display families of orbits for representative parameters.  Next we consider a more realistic predator-prey model that requires three independent parameters and supports more interesting phenomena, such as stable limit cycles.  In the process, we also introduce Mathematica functions for analyzing and plotting inequalities.

References:

J.D. Murray, Mathematical Biology, (Springer-Verlag, Berlin, 1989)

S. H. Strogatz, Nonlinear Dynamics and Chaos, (Perseus Books,  Reading MA, 1994)

Initialization

ClearAll["Global`*"];
Off[General::spell, General::spell1];

$DefaultFont={"Times", 12}; 
$TextStyle={FontFamily->"Times", FontSize->12,FontSlant->"Italic"};

Needs["Utilities`Notation`"];
Needs["Utilities`FilterOptions`"];
Needs["Graphics`"];
Needs["Algebra`InequalitySolve`"]

Original Lotka-Volterra model

The original Lotka-Volterra equations

d(X/d)τ == (A - B Y) 
        X
d(Y/d)τ == (C X - D) 
        Y

appear to involve 4 parameters, but by judicious scaling of X, Y, and τ we may be able to deduce simpler equations with fewer parameters.  Thus, let

X -> a x ,  
    Y -> b y,  τ -> c t

such that

d(x/d)t == (A  - b B y) 
        c x
d(y/d)t == (C a x - D) 
        c y

If we now choose

c  =  (A^(-1)),   b  =  A (B^(-1)),   a  =  
      D (C^(-1))

and define

α  =  D (A^(-1))

we obtain a simpler set of equations

d(x/d)t  ==  (1  -  y) 
        x
d(y/d)t  ==  α (x - 
            1) y

in terms of a single parameter, α, that governs the interaction between the two species.  Note that α > 0.

Although the fixed points are obviously {0, 0} and {1, 1}, we perform a derivation below to illustrate the steps.  First, we define the phase-space coordinates r = {x, y} and the corresponding velocity field v using the equations of motion.

r = {x, y};
v = {(1 - y) x, α (x - 1) y};

We then solve v[r] == 0 for the locations of the fixed points

fp = Solve[Thread[v == 0], r]

and construct the stability matrix S = (Subscript[((∂v/∂r)), fp]).

S = Outer[D, v, r] /. fp

Finally, we evaluate the eigenvalues of the stability matrix.

Eigenvalues /@ S

Thus, the origin has mixed stability: the prey population increases exponentially if there are no predators whereas the predator population decreases exponentially if there are no prey.  The fixed point at {1, 1} is neutral — the counterclockwise rotation about the fixed point represents stable oscillations in both populations.  Orbits close to the neutral fixed points are simply counterclockwise circles.  The same qualitative features are found for any positive α.

The following functions produce phase-space trajectories and plot them together with the velocity field.  Notice that we can supply a list of initial conditions, each a {(Subscript[x, 0]), (Subscript[y, 0])} pair, to PhaseAndField in order to plot a family of orbits on the same figure.  Also notice that we use FilterOptions to distribution the appropriate options to NDSolve, ParametricPlot, and PlotVectorField.  However, because some of the commonly used options for the two plotting functions are the same, some of the options for PlotVectorField were set to Automatic to avoid conflicts.  Some examples follow.

Clear[trajectory];
trajectory[α_, {x0_, y0_}, tmax_, 
        opts___Rule] := {x[t], 
          y[t]} /. NDSolve[{(x^′)[t] == (1 - y[t]) 
                  x[t], (y^′)[t] == α (x[t] - 1) 
                  y[t], x[0] == x0, y[0] == y0}, {x, y}, {t, 0, 
              tmax}, opts][[1]]

Clear[PhaseAndField];
\(PhaseAndField[α_, ic : {{_, _} .. }, xmax_, ymax_, tmin_, tmax_, 
        opts___Rule] := 
      With[{optParametricPlot\  = \ FilterOptions[ParametricPlot, opts], 
          optPlotField = FilterOptions[PlotVectorField, opts], 
          optND\  = 
            FilterOptions[NDSolve, 
              opts]}, \ 
DisplayTogether[
PlotVectorField[\ {\((1 - y)\) x, α\ \((x - 1)\) y}, {x, 0, xmax}, {y, 0, ymax}, 
            PlotPoints -> Automatic, 
            optPlotField], 
ParametricPlot[
            Evaluate[
              Table[trajectory[α, 
                  ic[[i]], tmax, 
                  optND], {i, Length[ic]}]], {t, tmin, tmax}, 
            optParametricPlot], Frame -> True, 
          FrameLabel -> {"prey", "predator"}, 
          PlotLabel -> "Lotka-Volterra model, α=" <> 
              ToString[α]]]\)

PhaseAndField[2, Table[{x0, 1}, {x0, 0.5, 4, 0.5}], 4, 6, 0, 
        10];

Plot[Evaluate[trajectory[2, {0.5, 0.5}, 200]], {t, 0, 10}, 
        PlotRange -> {Automatic, {0, 3}}, 
        PlotStyle -> {{}, Dashing[{0.02, 0.02}]}, Frame -> True, 
        FrameLabel -> {"t", "population density"}, 
        PlotLegend -> {"prey", "predator"}, 
        LegendPosition -> {-0.75, 0.4}, 
        LegendSize -> {0.7, 0.2}, 
        LegendShadow -> None];

We observe that the orbits are closed, corresponding to periodic oscillations for both populations.  These oscillations are approximately sinusoidal very close to the fixed point but become increasingly distorted as their amplitude increases.  Nevertheless, the orbits appear to remain closed for arbitrarily large amplitudes.  This is a new feature that is not predicted by our linear stability analysis, which applies in the immediate vicinity of the fixed point.  Why are the orbits closed?  It turns out that we can integrate these equation of motion for y[x] explicitly and thereby deduce a conserved quantity that is determined by the initial conditions and remains constant thereafter.  First, we eliminate the time dependence to obtain

d(y/d)x  ==  α ((x - 
              1)/(1 - 
              y)) (y/x)  ⟹   ((1 - 
                y)/y) dy  ==  -α ((1 - 
              x)/x) dx

Next we integrate both sides, such that

Log[y]  -  y  ==  -α (Log[x]   -  x)  -  
        T

where T is a constant.  Therefore, for specific choices of α and T we have a unique parametric dependence of y upon x; in other words, a closed orbit.  The constant of integration is determined by the initial conditions, such that

\(T\  = \ \(\(α\ x\_0\)\(\ \)\(-\)\(\ \)\(y\_0\)\(\ \)\(-\)\(\ \)\(Log[\(x\_0\^α\) y\_0]\)\(\ \)\)\)

is a constant of the motion.  Notice that T ≥ 1 + α and that the minimum value of T corresponds to the fixed point (Subscript[x, 0]) = (Subscript[y, 0]) = 1.  Furthermore, the greater the oscillation amplitude the larger the value of T.  The population parameter α defines a family of closed orbits while the conserved quantity T, as determined by the initial conditions, selects the appropriate member of that family.

Another feature of this system also merits discussion.  Notice that for large amplitudes the orbits tend to "bunch up" near the coordinate axes and in the lower left corner of phase space.  Thus, small fluctuations in environmental conditions could shift the system from one orbit to another nearby orbit whose amplitude is quite different in the region of phase space (here toward the upper right) where the orbits spread out again.  Therefore, small perturbations that occur at a vulnerable time can produce large effects upon the subsequent behavior of the system.  This is a common feature of conservative systems that is often not very realistic.

More realistic predator-prey model

More realistic predator-prey models should include saturation mechanisms for both species.  For example, we can describe the finite carrying capacity of the environment by assuming that the prey population would obey a logistic equation in the absence of predators.  Furthermore, we assume that there is maximum prey consumption rate per predator.  A simple equation with those features is

d(X/d)τ == (Subscript[R, X]) 
          X (1 - 
              (X/(Subscript[K, X])))  -  ((A X Y)/(X + B))

where (Subscript[R, X]) is a reproduction rate, (Subscript[K, X]) describes the environmental carrying capacity, A/B describes the predation rate for small prey populations, and B is a saturation parameter.  Similarly, we assume that the predator population is also described by a logistic equation but that the carrying capacity for predators is proportional to the prey density, such that

d(Y/d)τ == (Subscript[R, Y]) 
        Y (1 - (Y/((Subscript[K, Y]) X)))

As before, we would like to use dimensional analysis to reduce the number of free parameters.  Let

X -> α x ,  
    Y -> β y,  τ -> γ t
\(dX\/dτ == \(R\_X\) 
          X\ \((1 - 
              X\/K\_X)\)\  - \ \(\(A\ X\ Y\)\/\(X + 
                  B\)\)\ \ ⟹\ \ \(1\/\(γ\ x\)\) dx\/dt == \ 
      R\_X\ \((1\  - α\/K\_X\ x)\)\  - \ \(A\ β\ y\)\/\(α\ x\  + \ B\)\)
\(dY\/dτ == \(R\_Y\) 
        Y\ \((1 - 
              Y\/\(\(K\_Y\) X\))\)\ \ ⟹\ \ \(1\/\(γ\ y\)\) dy\/dt == \ \(\(R\_Y\)\(\ \)\((1 - \(β\/\(α\ K\_Y\)\) y\/x)\)\(\ \)\)\)

The predator equation takes its simplest form when we choose

\(β\/\(α\ K\_Y\) = 1\ , \ \ γ\ R\_Y = 
      c\ \ \ ⟹\ \ \(dy\/dt\)\  == \ 
        c\ y\ \((1 - y\/x)\)\)

Then the prey equation can be simplified by choosing

\(γ\  = R\_X\^\(-1\), \ \ α = 
      K\_X, \ \(A\ β\)\/K\_X = a, \ 
    B\/K\_X = 
      b\ \ \ ⟹\ \ \(dx\/dt\)\  == \ 
        x\ \((1 - x)\)\  - \ \(a\ x\ y\)\/\(x + b\)\)

Therefore, we obtain dimensionless equations of the form

d(x/d)t  ==  
      x (1 - x)  -  ((a x y)/(x + b))
d(y/d)t  ==  
      c y (1 - (y/x))

where

a  =  A (Subscript[K, Y]) ,   b  =  B (Subscript[K, X]) ,   c  =  
      ((Subscript[R, Y])/(Subscript[R, X]))

are three independent parameters that describe the dynamics of this system.  These parameters must be positive.

Next we evaluate the fixed points and the stability matrix for this system.

r = {x, y};
v = {x (1 - x)  -  ((a x y)/(x + b)), 
          c y (1 - (y/x))};

fp = Solve[Thread[v == 0], r]

S = Outer[D, v, r]

The second solution is not a viable fixed point with positive nonnegative populations.  The fixed point at {1, 0} is unstable whenever predators are present.

Eigenvalues[
      S /. fp[[1]]]

The nature of the nontrivial fixed point depends upon the parameter values in an interesting manner.  To analyze this behavior we observe that the secular equation for two-dimensional systems takes the form

Det[S - λ 𝟙] == 
      0  ⟹  (λ^2)  -  
        Tr[S] λ  +  Det[S]  ==  0

such that the eigenvalues

λ  =  (1/2) ( 
          Tr[S]  ±  (Sqrt[Tr[S]^2 - 4  Det[S]]) )

are expressed in terms of two simple properties of the stability matrix.  We assume that S is strictly real at the fixed point.  Therefore, the necessary and sufficient conditions for a stable fixed point require both Tr[S] < 0 and Det[S] > 0.  For this system these quantities become

trS = 
      FullSimplify[
        Tr[S /. fp[[3]]], {a > 0, b > 0, 
          c > 0}]

detS = 
      FullSimplify[
        Det[S /. fp[[3]]], {a > 0, 
          b > 0, c > 0}]

Although Mathematica has some difficulty analyzing this inequality for Det[S],

InequalitySolve[detS < 0, c] // 
      Simplify[#, {a > 0, b > 0}] &

one quickly sees that it cannot be satisfied unless c < 0, which is not a valid parameter range; therefore, the condition on Det[S] never leads to instability.  The condition on Tr[S], on the other hand, does yield an unstable region of parameter space.  This fixed point is stable for any a < (1/2), but for larger a and small b it is unstable when c is below the surface plotted below.

UnstableRegion = (Simplify[#, {a > 0, b > 0}] &) /@ 
        InequalitySolve[trS > 0, c]

InequalityPlot3D[UnstableRegion, {a, 0, 2}, {b, 0, 5}, {c, 0, 2}, 
        AxesLabel -> {"a", "b", "c"}, 
        ViewPoint -> {3.0,  2.0,  1.0}, 
        PlotLabel -> "unstable region: c below surface"];

The following functions evaluate the population history and produce figures that combine the vector field with orbits in phase space.

Clear[trajectory];
trajectory[{a_, b_, c_}, {x0_, y0_}, tmax_, 
        opts___Rule] := {x[t], 
          y[t]} /. NDSolve[{(x^′)[
                  t] == (1 - x[t]) 
                    x[t] - ((a x[t] y[t])/(x[t] + b)), (y^′)[
                  t] == c  y[t] (1 - y[t]/x[t]), 
              x[0] == x0, y[0] == y0}, {x, y}, {t, 0, tmax}, 
            opts][[1]]

Clear[PhaseAndField];
\(PhaseAndField[{a_, b_, c_}, ic : {{_, _} .. }, xmax_, ymax_, tmin_, 
        tmax_, opts___Rule] := 
      With[{optParametricPlot\  = \ FilterOptions[ParametricPlot, opts], 
          optPlotField = FilterOptions[PlotVectorField, opts], 
          optND\  = 
            FilterOptions[NDSolve, 
              opts]}, 
DisplayTogether[
PlotVectorField[\ {\((1 - x)\) x\  - \ \(a\ x\ y\)\/\(x + b\), 
              c\ \((1 - y\/x)\) y}, {x, 0.01, xmax}, {y, 0, ymax}, 
            PlotPoints -> Automatic, 
            optPlotField], 
ParametricPlot[
            Evaluate[
              Table[trajectory[{a, b, c}, 
                  ic[[i]], tmax, 
                  optND], {i, Length[ic]}]], {t, tmin, tmax}, 
            optParametricPlot], Frame -> True, 
          FrameLabel -> {"prey", "predator"}, 
          PlotLabel -> "predator-prey model, {a,b,c}=" <> 
              ToString[{a, b, c}]]]\)

Thus, we observe that all trajectories for a wide range of initial conditions rapidly spiral into a stable fixed.

ics = 
        Flatten[Table[{x0, y0}, {x0, 0.25, 1, 0.25}, {y0, 0.25, 1, 0.25}], 
          1];
PhaseAndField[{0.5, 0.1, 0.2}, ics, 1, 1, 0, 100];

Conversely, the figure below shows that trajectories for initial conditions near an unstable fixed point spiral out while those for more distant initial conditions spiral in toward a band that is suggestive of a common asymptotic orbit.

ics = 
        Flatten[Table[{x0, y0}, {x0, 0.1, 0.5, 0.1}, {y0, 0.1, 0.5, 0.1}], 
          1];
PhaseAndField[{1.5, 0.1, 0.2}, ics, 1, 0.5, 0, 100];

In fact, all of these orbits coincide for large times.  The initial conditions affect the phase on the asymptotic orbit but not its shape.

PhaseAndField[{1.5, 0.1, 0.2}, ics, 0.8, 0.35, 200, 400, 
        PlotPoints -> 50];

Finally, the figure below shows the time dependence of both population densities for a particular choice of initial conditions.  The greater amplitude of the oscillation of prey population and the approximately quarter-cycle relative phase between the two populations both make sense intuitively — the predator population peaks after the prey that it relies on for food, but increasing predator population drives down the prey population.  The relative scarcity of prey then drives down the predator population.  When the predator population is small, the prey population tends to increase.  For this system, the two populations eventually settle upon a periodic oscillation known as a limit cycle.

Plot[
        Evaluate[trajectory[{1.5, 0.1, 0.2}, {0.2, 0.2}, 200]], {t, 0, 80}, 
        PlotRange -> {Automatic, {0, 0.8}}, 
        PlotStyle -> {{}, Dashing[{0.02, 0.02}]}, Frame -> True, 
        FrameLabel -> {"t", "population density"}, 
        PlotLegend -> {"prey", "predator"}, 
        LegendPosition -> {-0.7, 0.35}, 
        LegendSize -> {0.7, 0.2}, 
        LegendShadow -> None];

Limit cycles require two conditions.  First, the fixed point must be unstable such that points inside the cycle are driven outward.  Second, there must be a closed surface in phase space on which the velocity field is everywhere inward.  When these conditions are satisfied, a limit cycle around the fixed point exists within that surface.  Thus, points inside the limit cycle are driven outward, while points outside are driven inward.  The limit cycle represents a stable periodic behavior — small perturbations decay and the system returns to the limit cycle, shifted only in phase. Limit cycles can be found in a wide range of phenomena, including population dynamics, chemical oscillators, mechanical or electrical systems, etc.

We have assumed that the model is autonomous — that there is no explicit time dependence.  More realistically, the model parameters may have small but rapid fluctuations or may gradually change with time due to external or environmental changes.  If the parameters change gradually without affecting the stability of the fixed point, we might expect the shape of the limit cycle to also change gradually (adiabatically) without affecting any qualitative features of the dynamics.  However, if the parameters of the system are near the transition between Re[λ] > 0 and Re[λ] < 0, fluctuations or parameter drift could have profound consequences: a stable system could suddenly begin oscillating or a periodic system could suddenly spiral toward a fixed point.