May 20, 2018, Sunday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Three Species Food Chain Modeling

From MathBio

Jump to: navigation, search

William Dowling, Toni Klopfenstein, Lucas Koepke


Songjuan Lv and Min Zhao examine a three-species food chain with a Holling-type functional response in the paper The Dynamic Complexity of a Three Species Food Chain Model [1]. By analytically examining and numerically simulating the system, the authors analyze the food chain’s chaotic, stable and periodic dynamics. The authors examine the four non-negative equilibrium points for the system, and perform a local stability analysis for each equilibrium point. The long-term behavior of the system is also examined using numerical analysis. It is also shown how the system can be pushed into chaos or into a steady state by varying certain parameters, such as the prey carrying capacity or the prey growth rate. An examination of the largest Lyapunov exponent is also included.


The three species food chain model is an extension of the general two species model given by Lotka and Volterra. Until recently, it was believed that only steady states and limiting cycles were exhibited in those types of systems. However, researchers have recently discovered the appearance of chaos and cycles in predator-prey species models with nonlinear responses, including Beddington-DeAngelis [2] and Holling-type responses.

This paper deals primarily with the Holling functional response. C.S. Holling[3] [4]discovered this response to describe the consumption rate of prey and the relationship with the prey population density. There are three different characteristic responses included in the Holling functional response family.[5]

In Type I responses, prey consumption grows linearly until reaching an equilibrium point at which consumption maintains a constant value. The consumption rate increases at the rate of predator attack rates. These are usually found in a herbiovoreoplant or invertebrate predator-prey food chains[6]

In Type II responses, prey consumption grows with an increase in prey density until a constant consumption level is reached, regardless of any other increases in prey density. This is the most commonly observed functional response in nature, and is modeled by Holling's Disk Equation, [7] shown below,

GeneralHolling.jpg where: HollingParameters.jpg

(This asymptotic response is more commonly seen in terms of N, as in APPM4390:Host-Parasitoid_Models#Type_II:_Asymptotic_in_N, and is also described in discrete time by the Holling Model - APPM4390:Host-Parasitoid_Models#Holling_Model)

Type III responses have an initial increase in prey consumption which then levels out to a constant consumption level.

The system examined in this paper involves a three-species trophic system, in which there are two competing predators (one intermediate and one top predator) and one prey species, and has a Holling Type-II functional response. Examples of this type of system include mouse-snake-owl, worm-robin-falcon, and vegetation-hare-lynx. [8]

Data taken from wild food chains is often incomplete and so it can be difficult to verify models. However in certain cases it is possible to reconstruct the states of all the species in a three species model. For more about food chain modeling see APPM4390:Robustness Analysis of an Observer-Based Controller in a Food Web

Mathematical Model

The general system of equations used to represent a three species food chain is:


given by Varriale and Gomes.[9] The parameters are defined as: Initialparameters.jpg

In this equation, the functions F1 and F2 are representations of the interactions between species, which represent a Holling-type II functional response, and are given by the equations below.


If the original system is modified to include this Holling-type functional response, the following model is achieved.


The new parameters are:


The model can then be transformed using:


and the dimensionless model is created:


with the new parameters:


Equilibria Analysis

Stability analysis of equilibria is performed on the non-dimensionalized form of the system given above.

There are four possible non-negative equilibrium points. They correspond to the following biological interpretations: (1) all species extinct; (2) prey population at carrying capacity in absence of both predators; (3) top predator absent and intermediate predator survives on its prey; (4) all species present.

The equilibrium points (1) and (2) corresponding to (0,0,0) and (k0,0,0) always exist. The Jacobian matrices of both are given below:

Example14.jpg Example15.jpg

The eigenvalues of both are given by the diagonal entries in the corresponding Jacobian matrix. Both have one non-negative and two negative eigenvalues, from which we determine that both equilibria are saddle points. (Gakkhar et al similarly found saddle point equilibria when looking at delayed prey-predator models: APPM4390:AliHussain#Result)

The equilibrium point (x*,y*,0) corresponding to (3) where x* and y* are defined as


exists only when the following conditions are met:


In order to study the stability of (3), we can study the following two dimensional subsystem for which (3) is an equilibrium point:


The Jacobian of this system at (x*,y*) is given by:

Example17a.jpg where


This equilibrium is asymptotically stable if both coefficients of the corresponding characteristic equation are greater than zero. The coefficients for the characteristic equation


are given by


The conditions given for the existence of an equilibrium at (x*,y*,0) force the second coefficient to be greater than zero. Below are two graphs showing the change in the first coefficient with respect to r0 and k0:

Graph3r0.jpg Graph3k0.jpg

Note that for the range of values tested, the value of the first coefficient never surpasses zero except for k0 < 1.3 (roughly). As a result, we can determine that the equilibrium point is not asymptotically stable for all values tested except for k0 < 1.3. However, the condition 0 < x* < k0 < k1 must still hold for the equilibrium point to exist (for the parameter used here, x* = 1/3).

The fourth equilibrium exists for positive x, y, and z values if and only if there exists a positive solution to the following system of equations (derived directly from the non-dimensionalized system of differential equations being studied):


The solution is given by

Example10.jpg, Example11.jpg, and Example12b.jpg


Example12a.jpg and Example12c.jpg

The Jacobian at (4) is given by




The coefficients of the corresponding characteristic equation


are given by


The equilibrium point (4) is locally stable if the following condition is met:


Below, the quantities Sigma1.jpg and Sigma2.jpgare graphed for a range of r0 and k0 values (the median possible x value is used):

Graph4r0.jpg Graph4k0.jpg

Note that all must be greater than zero in order for the equilibrium point to be asymptotically stable. This does not occur for any of the r0 or k0 values tested, thus we can determine that the equilibrium (4) is not asymptotically stable for any of these values.

Numerical Analysis

To study the behavior of the non-dimensionalized system, we will examine how varying one of the input parameters, while keeping the others fixed, affects the system.

As an example, we look at the parameter r0, which is the growth rate of species x, and how it determines the characteristics of the top predator z. To do this, we will use a Poincare map of the local maxima for the top predator z. The fixed parameter values are: PoincareParameters.jpg with initial condition x0 = 3, y0 = 0.1 and z0 = 0.2

Poincare map

The Poincare map is a useful tool in the analysis of dynamical systems, and is defined as follows. Given an n-dimensional system, let S be an n-1 dimensional surface such that all flow through S is entirely perpendicular to it. That is, there are no trajectories, starting on S, with flow parallel to S. Then the Poincare map, call it P, is defined as a mapping from S to S, such that

xk+1 = P (xk)


Then a closed orbit of the original system will be a fixed point of P, such that after some time t, P (x*) = x*. In practice, it is nearly impossible to explicitly determine this mapping P, however using numerical techniques the Poincare map provides a useful framework to determine when cycles exist in a dynamical system. Specifically, the continuous solution to the original equation will have to be discretized in some way to form a sequence a values that can be analyzed for points that fit the definition for a fixed point for P. A natural way to do this is to look at the sequence of values corresponding to the local maxima, and any values that repeat will be fixed points of P.

Numerical results

To use this technique on the non-dimensionalized system given above, we will let r0, which is the growth rate of the prey x, vary between 1.5 and 9, and then analyze the changes in behavior produced in the population of the top predator z. A sequence of the local maxima of the time series for z is generated, and then any repeated values in this sequence show a cycle exists in the original system.

To begin with, for r0 = 1.5, the plot of z with respect to time looks like


and highlighting the maxima with red dots looks like


Starting from the left side of the picture, we let x1 be the first maxima, and then x2 = P(x1) be the second, and so forth. The sequence of values generated for the local maxima is then

{1.416, 1.076, 1.207, 1.199, 1.201, 1.201, 1.2, 1.201, 1.201, 1.201, 1.201, 1.201}

And since 1.201 is repeated, and is in fact the only repeated value, we now have a single fixed point for P, which corresponds to a period-1 cycle in the original system.

Similarly, for r0 = 3, the plot looks like


and highlighting the maxima with red dots looks like


Then the sequence of values generated for the local maxima is

{1.909, 0.957, 2.091, 0.778, 2.082, 0.785, 2.081, 0.786, 2.081, 0.787, 2.081, 0.787, 2.081, 0.787, 2.081}

which has two fixed points for P, which correspond to a period-2 cycle in the original system.

To illustrate the dynamics better, we let r0 vary between 1.5 and 9, and plot only the fixed points of the Poincare map. This will show where bifurcations in the system occur, such as a period-1 cycle splitting into a period-2 cycle, this period-2 cycle bifurcating again into a period-4 cycle, etc. The resulting image, with z on the vertical axis and r0 on the horizontal axis, is


As r0 is further increased, we see period doubling, first to a period-4 cycle, and then chaos until r0 = 6.5 where a period-3 cycle emerges, only to undergo period doubling again, leading to chaos.


As shown by the stability analysis, we determined that separately varying two of the parameters, r0 and k0, within the set limits would never allow an asymptotically stable solution. This was first shown using the Routh-Hurwitz criteria for the coefficients in the characteristic equation, and since one of them was always less than zero, no asymptotically stable solution could exist. This was verified with the numerical results for r0 between 1.5 and 9, since there was always at least one fixed point for the Poincare map, and hence there was always at least one cycle in the original system. This periodic nature of the solutions to our model for three interacting species fits our intuition, since we would not expect to have an exactly constant number of any one species in the real world.

Since we were able to reproduce and verify the results found in the paper we reviewed, using Matlab and Mathematica, we conclude that a justifiable and accurate model was used. The methods for both stability analysis and numerical modeling are techniques commonly applied to dynamical systems, and the results from both of these methods showed the same conclusion. The occurrence of bifurcations and chaos are also possible with a three species model, and both of these cases were observed using numerical modeling. As a result, studying a three species food chain yields complex dynamics and interesting behavior, lending itself well to mathematical insight, while still staying grounded in reality.

Project categorization


This model uses a system of ordinary differential equations to model the population dynamics of three interacting species. Stability analysis of the found equilibrium points is performed using Jacobian matrices and eigenvalues.

Model type

The system of ODEs are modeling the change in each of the different species populations with respect to time.

Biological system

The biological system studied is the interaction of three competing species.

Related Topics

Another example of useful population modeling techniques would be the use of Trophic pyramids which don't count individuals, but rather net biomass or production rates.

Recent Citation

The above article was recently cited in Min Zhao and Limin Zhang’s 2009 paper titled “Permanence and chaos in a host–parasitoid model with prolonged diapause for the host”. The paper discusses the effects of a prolonged dormant period by the parasitoid, and how it alters the behavior of the host parasitoid relationship. The relationship model is described by a pair of coupled difference equations.

The permanence of the system is investigated through bifurcation analysis and using the largest Lyapunov exponents while altering the diapause rate d. It was found that when 0.281 < d < 0.4879, the Lyapunov exponents are negative, which means that there exist a stable coexistence state between the host and parasitoid.

The article was cited when looking at the complexities of the observed system. The parasitoid host model has similar complex dynamic patterns found in the bifurcation analysis as the three – species food-chain model.


  1. The dynamic complexity of a three species food chain model
  2. J. R. Beddington, Mutual Interference Between Parasites or Predators and its Effect on Searching Efficiency, Journal of Animal Ecology, Vol. 44, No. 1 (Feb., 1975), pp. 331-340
  3. Holling, C. S. 1959. Some characteristics of simple types of predation and parasitism. Canadian Entomologist 91:385-398
  4. Holling, C.S.: The components of predation as revealed by a study of small mammal predation of the European pine sawfly. Can. Entomologist 91, 293–320 (1959)
  5. Functional Responses
  6. Type I Functional Response
  7. Type II Functional Response: Holling's Disk Equation
  8. A Lotka-Volterra Three-Species Food Chain
  9. A study of a three species food chain