Table of Contents
ISRN Applied Mathematics
Volume 2012, Article ID 391547, 11 pages
Research Article

A Mathematical Model of Three-Species Interactions in an Aquatic Habitat

Department of Mathematics, University of Jos, PMB 2084, Jos, Nigeria

Received 16 November 2011; Accepted 13 December 2011

Academic Editors: M. Mei and X. Meng

Copyright © 2012 J. N. Ndam et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


A mathematical model for three-species interactions in a food chain, with the assumption that the interacting species are mobile, has been constructed using a combination of Holling’s type III and the BD functional responses. Conditions for the onset of diffusive instability were determined. The results indicate the possibility of a stable coexistence of the three interacting species in form of stable oscillations under the reflecting boundary conditions. Habitat segregation also occurs under these conditions. However, under the absorbing boundary conditions, the species experience damped oscillations leading to their extinction. The effects of cross-diffusion of the intermediate and the toppredator were also examined.

1. Introduction

A food chain is a common natural occurrence in an ecosystem due to interactions among species in their habitats. A food chain could occur either among terrestrial or aquatic species. The relationships among species could be those of predator-prey, mutualism, or competition for resources. In this paper, we consider predator-prey interactions in a food chain. From the Lotka-Volterra food chain model, a lot has been done on predator-prey interactions in food chains using different functional responses such as the Holling’s types (I–IV) and the Beddington-DeAngelis (BD) functional response and many others [13]. The Holling’s type II functional response is the most commonly used in mathematical ecology [4]. However, this type of functional response may not be the best for models involving interference among predators, as this is best modelled by the BD functional response which accommodates interference among predators [4, 5]. However, the Holing’s type III is more suitable for models which assume resource refuge effects, that is, scarce resources become better protected from exploitation. Kassem and Ndam [6] have exploited some of these special properties of the functional responses to construct a model for the predator-prey interactions in a food chain. We intend to build on this in the current model by assuming that the three species are mobile, which makes it more physically realistic, since species move about and also try to escape from their enemies in their natural habitats. Apart from exploiting those properties as explained by Kassem and Ndam [6], self-diffusion as well as cross-diffusion of the species will be examined, as can be seen in Okubo and Levin [7]. With these assumptions, we hope to capture most of the physical interactions of species in an aquatic or terrestrial habitat. The Holling’s type III functional response is given by𝑓(𝑃)=𝑐𝑃2𝑎2+𝑃2,(1.1) where 𝑎 is the half-saturation constant and 𝑐 is the maximal consumption rate of the predator and 𝑃 is the prey density, while the BD functional response is given by𝑓(𝑃,𝑄)=𝛼𝑃,𝜆+𝑃+𝛾𝑄(1.2) where 𝑄 is the intermediate predator density, while 𝛼, 𝜆, and 𝛾 are positive constants, representing the maximum consumption rate, the saturation constant, and the predator interference parameter, respectively [8]. The remaining parts of this paper will be organised as follows: the formulation of the model will be considered in Section 2, numerical simulations will be considered in Section 3, while Section 4 will be dedicated to the stability analysis of the system. Finally, some conclusions and future directions will be mentioned.

2. Mathematical Formulation

A model for predator-prey interactions in a three-trophic level food chain with diffusion of the three species is constructed. This model consists of a primary producer 𝑃, called the prey, an intermediate predator 𝑄, and a top predator 𝑅, a term adopted from [2]. We assume that the prey 𝑃 has a logistic growth in the absence of the intermediate predator. It is also assumed that the intermediate predator 𝑄 has a BD-functional response, while the top predator 𝑅, has a Holling’s type III functional response and that all the three are mobile. Hence we obtain the following model equations:𝜕𝑃𝜕𝑡=𝐷1𝜕2𝑃𝜕𝑥2𝑃+𝑟𝑃1𝐾𝑎1𝑃𝑄,𝜆+𝑃+𝛾𝑄(2.1)𝜕𝑄𝜕𝑡=𝐷2𝜕2𝑄𝜕𝑥2+𝑏1𝑃𝑄𝑎𝜆+𝑃+𝛾𝑄2𝑄2𝑅𝛽2+𝑄2𝜇1𝑄,(2.2)𝜕𝑅𝜕𝑡=𝐷3𝜕2𝑅𝜕𝑥2+𝑏2𝑄2𝑅𝛽2+𝑄2𝜇2+𝑞𝑅,(2.3)𝑃(𝑥,0)>0,𝑄(𝑥,0)>0,𝑅(𝑥,0)>0(2.4) subject to the reflecting boundary conditions𝐷1𝜕𝑃(0,𝑡)𝜕𝑥=𝐷1𝜕𝑃(𝐿,𝑡)𝜕𝑥=0,𝐷2𝜕𝑄(0,𝑡)𝜕𝑥=𝐷2𝜕𝑄(𝐿,𝑡)𝜕𝑥=0,𝐷3𝜕𝑅(0,𝑡)𝜕𝑥=𝐷3𝜕𝑅(𝐿,𝑡)𝜕𝑥=0(2.5) or the absorbing boundary conditions𝑃(0,𝑡)=𝑃(𝐿,𝑡)=0,𝑄(0,𝑡)=𝑄(𝐿,𝑡)=0,𝑅(0,𝑡)=𝑅(𝐿,𝑡)=0,(2.6) where 𝑟>0 is the intrinsic growth rate of the prey and 𝐾>0 is its carrying capacity and 𝐷1, 𝐷2, and 𝐷3 are the constant diffusion coefficients of the prey, the intermediate predator, and the top-predator, respectively; 𝜇1>0, 𝜇2>0, and 𝑞>0 are the natural mortality rates of the intermediate and top-predators, and the combined harvesting effort of the top-predator respectively. The other positive constants 𝑎𝑖, 𝑏𝑗, 𝛽, 𝜆, and 𝛾 have the same meanings as defined in (1.1) and (1.2). The reflecting boundary conditions (2.5) imply that there are no population fluxes across the boundaries of the habitat. The boundary conditions (2.6), on the other hand mean that the environment always contains the equilibrium population [7].

Equations (2.1)–(2.3) can be scaled as follows:𝑡𝑡=𝑟,𝑃=𝐾𝑃,𝑄=𝐾𝑄,𝑅=𝐾𝑅,𝑥=𝐿𝑥,𝑎1=𝑟𝑎1,𝑎2=𝑟𝑎2,𝑏1=𝑟𝑏1,𝑏2=𝑟𝑏2,𝜇1=𝑟𝜇1,𝜇2=𝑟𝜇2,𝑞=𝑟𝑞,(2.7) where the asteriated variables are nondimensional variables. We obtain after dropping asterisks the scaled system𝜕𝑃=𝜕𝜕𝑡2𝑃𝜕𝑥2𝑎+𝑃(1𝑃)1𝑃𝑄𝑚1,+𝑃+𝛾𝑄(2.8)𝜕𝑄𝜕𝜕𝑡=𝜌2𝑄𝜕𝑥2+𝑏1𝑃𝑄𝑚1𝑎+𝑃+𝛾𝑄2𝑄2𝑅𝑚22+𝑄2𝜇1𝑄,(2.9)𝜕𝑅𝜕𝜕𝑡=𝜎2𝑅𝜕𝑥2+𝑏2𝑄2𝑅𝑚22+𝑄2𝜇2+𝑞𝑅,(2.10)𝑃(𝑥,0)>0,𝑄(𝑥,0)>0,𝑅(𝑥,0)>0,(2.11)𝜕𝑃(0,𝑡)=𝜕𝑥𝜕𝑃(𝐿,𝑡)𝜕𝑥=0,𝜌𝜕𝑄(0,𝑡)𝜕𝑥=𝜌𝜕𝑄(𝐿,𝑡)𝜕𝑥=0,𝜎𝜕𝑅(0,𝑡)𝜕𝑥=𝜎𝜕𝑅(𝐿,𝑡)𝜕𝑥=0,(2.12) where we have taken 𝐿=𝐷1/𝑟, 𝜌=𝐷3/𝐷1,𝜎=𝐷2/𝐷1, 𝑚1=𝜆/𝐾, and 𝑚2=𝛽/𝐾.

2.1. The Effects of Cross-Diffusion

It is natural for the preys to run away in the presence of their enemies, the top-predators. This phenomenon is called cross-diffusion. We introduce in this model cross-diffusion terms with constant coefficients 𝐷12 and 𝐷21. The scaled model equations then take the form𝜕𝑃=𝜕𝜕𝑡2𝑃𝜕𝑥2𝑎+𝑃(1𝑃)1𝑃𝑄𝑚1,+𝑃+𝛾𝑄(2.13)𝜕𝑄𝜕𝜕𝑡=𝜌2𝑄𝜕𝑥2+𝜕2𝑅𝜕𝑥2+𝑏1𝑃𝑄𝑚1𝑎+𝑃+𝛾𝑄2𝑄2𝑅𝑚22+𝑄2𝜇1𝑄,(2.14)𝜕𝑅𝜕𝜕𝑡=𝜎2𝑅𝜕𝑥2𝜕+𝛿2𝑄𝜕𝑥2+𝑏2𝑄2𝑅𝑚22+𝑄2𝜇2+𝑞𝑅,(2.15) where 𝛿=𝐷21/𝐷12, subject to the boundary conditions (2.6) and (2.12).

3. Numerical Simulation

Numerical simulations showing the spatial diffusions and the time-course solution are shown in the figures below for the parameter values𝑎1=1.2,𝑎2=0.25,𝑏1=0.5,𝑏2=0.75,𝑚1=0.2,𝑚2𝜇=0.1,1=𝜇2𝑎=0.1,𝛾=0.25,𝜌=1.25,𝜎=1.05,𝑞=0.35,1=2.25,𝑎2=0.15,𝑏1=0.75,𝑏2=1.05,𝑚1=0.2,𝑚2𝜇=0.1,1=𝜇2=0.1,𝛾=0.5,𝜌=1.10,𝜎=1.5,𝑞=0.5,(3.1) for the two sets of boundary conditions (2.6) and (2.12). The results indicate that while the absorbing boundary conditions lead to diffusion of the three interacting species, the population density distributions depict a scenario known as habitat segregation under the reflecting boundary conditions as seen in Figures 1 and 2, respectively. On the other hand, the time-course solution under the absorbing boundary conditions produce damped oscillations leading to the gradual extinction of the interacting species (Figure 3). However, imposing the reflecting boundary conditions leads to a coexistence of the species in form of stable oscillations (Figure 4). Moreover, the long-time behaviour of the interacting species indicates the appearance of patchiness in such a way that the preys tend to separate from their predators (Figure 5). This result confirms the fact that habitat segregation serves as a stabilising factor for possible coexistence of the interacting species, as noted by Okubo and Levin [7].

Figure 1: Diffusion of the three species for the absorbing boundary conditions.
Figure 2: Habitat segregation due to the reflecting boundary conditions.
Figure 3: Extinction of the three species with the absorbing boundary conditions.
Figure 4: Stable oscillations of the interacting species with time.
Figure 5: Development of patchiness as a long time scenario.

Figures 6 and 7 show the effects of cross-diffusion on the spatial density distribution of the interacting species. While Figure 6 depicts the influence of both self- and cross-diffusion on the spatial distribution of the species where the intermediate predators diffuse faster than the others, Figure 7 shows the population distribution due to cross-diffusion in the absence of self-diffusion, in which case the top-predator becomes least mobile. This could be due to the habit segregation which occurs under the circumstance. The scenarios described above clearly explain what actually happens among interacting species in an ecosystem. In the presence of self- and cross-diffusion, both species are mobile, and in addition, the intermediate predators try to escape from the top-predators, hence they are most mobile.

Figure 6: Cross and self-diffusion of the interacting species.
Figure 7: Spatial distribution of population due to cross-diffusion of the predators.

4. Stability of the Uniform Steady State

In order to examine the stability, the system (2.13) can be recast in the form𝜕𝑃=𝜕𝜕𝑡2𝑃𝜕𝑥2+𝐹1(𝑃,𝑄),𝜕𝑄𝜕𝜕𝑡=𝜌2𝑄𝜕𝑥2+𝜕2𝑅𝜕𝑥2+𝐹2(𝑃,𝑄,𝑅),𝜕𝑅𝜕𝜕𝑡=𝜎2𝑅𝜕𝑥2𝜕+𝛿2𝑄𝜕𝑥2+𝐹3(𝑄,𝑅),(4.1) where𝐹1𝑎(𝑃,𝑄)=𝑃(1𝑃)1𝑃𝑄𝑚1,𝐹+𝑃+𝛾𝑄2𝑏(𝑃,𝑄,𝑅)=1𝑃𝑄𝑚1𝑎+𝑃+𝛾𝑄2𝑄2𝑅𝑚22+𝑄2𝜇1𝐹𝑄,3𝑏(𝑄,𝑅)=2𝑄2𝑅𝑚22+𝑄2𝜇2+𝑞𝑅.(4.2) Suppose that there exists a spatially uniform steady state (𝑃𝑐,𝑄𝑐,𝑅𝑐), such that𝐹𝑗𝑃𝑐,𝑄𝑐,𝑅𝑐=0,𝑗=1,2,3.(4.3) To determine the stability of the uniform steady state, we perturb the population densities thus𝑃(𝑥,𝑡)=𝑃𝑐+𝑃𝑄(𝑥,𝑡),(𝑥,𝑡)=𝑄𝑐+𝑄(𝑥,𝑡),𝑅(𝑥,𝑡)=𝑅𝑐+𝑅(𝑥,𝑡).(4.4) Hence the linearised system becomes𝜕𝑃=𝜕𝜕𝑡2𝑃𝜕𝑥2+𝑎11𝑃+𝑎12𝑄+𝑎13𝑅,𝜕𝑄𝜕𝜕𝑡=𝜌2𝑄𝜕𝑥2+𝜕2𝑅𝜕𝑥2+𝑎21𝑃+𝑎22𝑄+𝑎23𝑅,𝜕𝑅𝜕𝜕𝑡=𝜎2𝑅𝜕𝑥2𝜕+𝛿2𝑄𝜕𝑥2+𝑎31𝑃+𝑎32𝑄+𝑎33𝑅,(4.5) where𝑎𝑖𝑗=𝜕𝐹1𝜕𝑃𝜕𝐹1𝜕𝑄𝜕𝐹1𝜕𝑅𝜕𝐹2𝜕𝑃𝜕𝐹2𝜕𝑄𝜕𝐹2𝜕𝑅𝜕𝐹3𝜕𝑃𝜕𝐹3𝜕𝑄𝜕𝐹3𝜕𝑅(4.6) is the Jacobian matrix for the system evaluated at the uniform equilibrium point (𝑃𝑐,𝑄𝑐,𝑅𝑐). In order to examine the linear stability of the system (4.5), we use the normal modes approach by assuming solutions of the form 𝑃𝑒(𝜆𝑡+𝑖𝑘𝑥),𝑄𝑒(𝜆𝑡+𝑖𝑘𝑥), and 𝑅𝑒(𝜆𝑡+𝑖𝑘𝑥), where 𝜆 and 𝑘 are constants. Thus the eigenvalues satisfy the equation|||||||||𝜆+𝑘2𝑎11𝑎120𝑎21𝜆+(1+𝜌)𝑘2𝑎22𝑎230𝑎32𝜆+(𝜌+𝛿)𝑘2𝑎33|||||||||=0,(4.7) and one obtains the characteristic polynomial𝜆3+𝑑11+𝑑22+𝑑33𝜆2+𝑑11𝑑22+𝑑11𝑑33+𝑑22𝑑33𝑎12𝑎21𝑎23𝑎32𝜆+𝑑11𝑑22𝑑33𝑎23𝑎32𝑑11𝑎12𝑎21𝑑33=0,(4.8) where 𝑑11=𝑘2𝑎11,𝑑22=(1+𝜌)𝑘2𝑎22 and 𝑑33=(𝜌+𝛿)𝑘2𝑎33.

For stability, the following conditions should be satisfied: 𝑑11+𝑑22+𝑑33𝑑>0,(4.9)11𝑑22𝑑33𝑎23𝑎32𝑑11𝑎12𝑎21𝑑33𝑑>0,(4.10)11+𝑑22+𝑑33𝑑11𝑑22+𝑑11𝑑33+𝑑22𝑑33𝑎12𝑎21𝑎23𝑎32>𝑑11𝑑22𝑑33𝑎23𝑎32𝑑11𝑎12𝑎21𝑑33.(4.11)

If diffusion is neglected, that is, 𝑘=0,𝑑11+𝑑22+𝑑33>0𝑎11+𝑎22+𝑎33<0.(4.12)

Diffusive instability sets in if any of conditions (4.9)–(4.11) is violated. Condition (4.9) is invariably satisfied in the presence of diffusion. Thus, only conditions (4.10) and (4.11) will be examined. For instability in the system due to self-diffusion (𝛿=0), inequalities (4.10) and (4.11) are reversed, that is𝐻𝑘2=𝐴1𝑘6𝐴2𝑘4+𝐴3𝑘2+𝐴4𝐻<0,(4.13)𝑘2=𝐵1𝑘6𝐵2𝑘4+𝐵3𝑘2+𝐵4<0,(4.14)

where 𝐴1=𝜌+𝜌2,𝐴2=𝑎11𝜌2+(𝑎11+𝑎22+𝑎33)𝜌+𝑎33,𝐴3=(𝑎11𝑎22+𝑎11𝑎33𝑎12𝑎21)𝜌+𝑎11(𝑎22+𝑎33) and  𝐴4=𝑎11𝑎23𝑎32+𝑎12𝑎21𝑎33𝑎11𝑎22𝑎33;  while  𝐵1=(2𝜌+1)(𝜌+2)(𝜌+1),𝐵2=(2𝜌+1)(2𝜌+3)𝑎11+3(𝜌+1)𝑎222+(𝜌+2)(3𝜌+2)𝑎33,  𝐵3=(2𝜌+1)(𝑎211𝑎23𝑎32)+4(𝜌+1)(𝑎11𝑎22+𝑎11𝑎33+𝑎22𝑎33)+(𝜌+2)(𝑎233𝑎12𝑎21)+(𝜌+1)𝑎222, and  𝐵4=(𝑎11+𝑎22)(𝑎12𝑎21𝑎233)+(𝑎22+𝑎33)(𝑎23𝑎32𝑎211)(𝑎11+𝑎33)𝑎2222𝑎11𝑎22𝑎33.

The minimum value of 𝐻(𝑘2) occurs at 𝑘2=𝑘2min=(𝐴2+(𝐴223𝐴1𝐴3))/(3𝐴1).

Hence for instability to set in, 𝐻𝑘2min=𝐴2+𝐴223𝐴1𝐴3327𝐴21+𝐴3𝐴2+𝐴223𝐴1𝐴33𝐴1+𝐴4<𝐴2𝐴2+𝐴223𝐴1𝐴329𝐴21.(4.15)

Similar result is obtained for 𝐻(𝑘2min) in terms of the 𝐵𝑖,𝑖=1,2,3,4 in (4.14).

5. Conclusions

A mathematical model for three-species interactions in a food chain is constructed using Holling’s type III and the BD functional responses. Two sets of initial and boundary conditions, the absorbing and the reflecting boundary conditions, are imposed. The species are also assumed to undergo self- and cross-diffusion with constant coefficients of diffusion. Numerical simulation results show that the absorbing boundary conditions allow mobility of the species throughout the spatial domain. The top-predators tend to diffuse faster than the intermediate predators, even when the coefficient of diffusion of the latter is higher (Figure 1). This situation could be explained by the fact that the intermediate predators are protected at low densities by the use of the Holling’s type III functional response. However, under these boundary conditions, the species are driven into extinction as depicted in Figure 3.

With the reflecting boundary conditions, the habitat segregates and the species coexist in form of stable oscillations as can be seen from Figures 2 and 4, respectively. The long-time scenario indicates the development of patchiness in which the preys tend to separate from their predators (Figure 5). When both self- and cross-diffusions are present, the population density distributions indicate that the intermediate predator diffuse faster than the others as depicted in Figure 6. However, in the absence of self-diffusion, the top-predator becomes the least mobile (see Figure 7). The change in scenario is probably due to the habitat segregation. This model, with the reflecting boundary conditions, can be applied to fisheries, where the economically more viable species is harvested periodically without driving them into extinction.

Conditions for the onset of diffusive instability are also determined using linear stability theory. Nonlinear stability analysis becomes cumbersome and it is left for future investigation.


  1. S. B. Hsu, Y.-S. Li, and P. Waltman, “Competition in the presence of a lethal external inhibitor,” Mathematical Biosciences, vol. 167, no. 2, pp. 177–199, 2000. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  2. C.-H. Chiu and S.-B. Hsu, “Extinction of top-predator in a three-level food-chain model,” Journal of Mathematical Biology, vol. 37, no. 4, pp. 372–380, 1998. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  3. A. Basset, D. L. DeAngelis, and J. E. Diffendorfer, “The effect of functional response on stability of a grazer population on a landscape,” Ecological Modelling, vol. 101, no. 2-3, pp. 153–162, 1997. View at Publisher · View at Google Scholar
  4. G. T. Skalski and J. F. Gilliam, “Functional responses with predator interference: viable alternatives to the Holling type II model,” Ecology, vol. 82, no. 11, pp. 3083–3092, 2001. View at Google Scholar
  5. R. S. Cantrell and C. Cosner, “On the dynamics of predator-prey models with the Beddington-DeAngelis functional response,” Journal of Mathematical Analysis and Applications, vol. 257, no. 1, pp. 206–222, 2001. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  6. T. G. Kassem and J. N. Ndam, “On the dynamics of predator-prey interactions in a three-trophic level food chain,” Far East Journal of Applied Mathematics, vol. 34, no. 1, pp. 95–110, 2009. View at Google Scholar · View at Zentralblatt MATH
  7. A. Okubo and S. A. Levin, Diffusion and Ecological Problems: Modern Perspectives, Springer, New York, NY, USA, 2001.
  8. S.-B. Hsu, T.-W. Hwang, and Y. Kuang, “A ratio-dependent food chain model and its applications to biological control,” Mathematical Biosciences, vol. 181, no. 1, pp. 55–83, 2003. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet