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

Abstract

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 [1–3]. 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𝑓(𝑃)=𝑐𝑃²𝑎²+𝑃²,(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.1)𝜕𝑡=𝐷²𝜕²𝑄𝜕𝑥²+𝑏¹𝑃𝑄𝜆+𝑃+𝛾𝑄−𝑎²𝑄²𝑅𝛽²+𝑄²𝜕𝑅−𝜇₁𝑄,(2.2)𝜕𝑡=𝐷³𝜕²𝑅𝜕𝑥²+𝑏²𝑄²𝑅𝛽²+𝑄²𝑃(𝑥,0)>0,𝑄(𝑥,0)>0−,𝜇₂𝑅+(𝑥𝑞,0𝑅),>0((22..34)) subject to the reflecting boundary conditions𝐷₁𝜕𝑃(0,𝑡)𝜕𝑥=𝐷¹𝜕𝑃(𝐿,𝑡)𝜕𝑥=0,𝐷²𝜕𝑄(0,𝑡)𝜕𝑥=𝐷²𝜕𝑄(𝐿,𝑡)𝜕𝑥=0,𝐷³𝜕𝑅(0,𝑡)𝜕𝑥=𝐷³𝜕𝑅(𝐿,𝑡)𝜕𝑥=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 𝐷₁, 𝐷₂, and 𝐷₃ are the constant diffusion coefficients of the prey, the intermediate predator, and the top-predator, respectively; 𝜇₁>0, 𝜇₂>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,𝑎₂=𝑟𝑎^∗2,𝑏₁=𝑟𝑏^∗1,𝑏₂=𝑟𝑏^∗2,𝜇₁=𝑟𝜇^∗1,𝜇₂=𝑟𝜇^∗2,𝑞=𝑟𝑞^∗,(2.7) where the asteriated variables are nondimensional variables. We obtain after dropping asterisks the scaled system𝜕𝑃𝜕𝑡=𝜕²𝑃𝜕𝑥²+𝑃(1−𝑃)−𝑎¹𝑃𝑄𝜕𝑄𝑚₁+𝑃+𝛾𝑄,(2.8)𝜕𝑡=𝜌𝜕²𝑄𝜕𝑥²+𝑏¹𝑃𝑄𝑚₁+𝑃+𝛾𝑄−𝑎²𝑄²𝑅𝑚²²+𝑄²𝜕𝑅−𝜇₁𝑄,(2.9)𝜕𝑡=𝜎𝜕²𝑅𝜕𝑥²+𝑏²𝑄²𝑅𝑚²²+𝑄²𝜕𝑃(0,𝑡)𝑃(𝑥,0)>0,𝑄(𝑥,0)>0−,𝜇₂𝑅+(𝑥𝑞,0𝑅),>0,((22..1101))𝜕𝑥=𝜕𝑃𝜕(𝑥𝐿,𝑡)=0,𝜌𝜕𝑄𝜕(𝑥0,𝑡)=𝜌𝜕𝑄𝜕(𝑥𝐿,𝑡)=0,𝜎𝜕𝑅𝜕(𝑥0,𝑡)=𝜎𝜕𝑅𝜕(𝑥𝐿,𝑡)=0,(2.12) where we have taken 𝐿=√𝐷₁/𝑟, 𝜌=𝐷₃/𝐷₁,𝜎=𝐷₂/𝐷₁, 𝑚₁=𝜆/𝐾, and 𝑚₂=𝛽/𝐾.

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 𝐷₁₂ and 𝐷₂₁. The scaled model equations then take the form𝜕𝑃𝜕𝑡=𝜕²𝑃𝜕𝑥²+𝑃(1−𝑃)−𝑎¹𝑃𝑄𝜕𝑄𝑚₁+𝑃+𝛾𝑄,(2.13)𝜕𝑡=𝜌𝜕²𝑄𝜕𝑥²+𝜕²𝑅𝜕𝑥²+𝑏¹𝑃𝑄𝑚₁+𝑃+𝛾𝑄−𝑎²𝑄²𝑅𝑚²²+𝑄²𝜕𝑅−𝜇₁𝑄,(2.14)𝜕𝑡=𝜎𝜕²𝑅𝜕𝑥²+𝛿𝜕²𝑄𝜕𝑥²+𝑏²𝑄²𝑅𝑚²²+𝑄²−𝜇₂+𝑞𝑅,(2.15) where 𝛿=𝐷₂₁/𝐷₁₂, 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.2,𝑎₂=0.25,𝑏₁=0.5,𝑏₂=0.75,𝑚₁=0.2,𝑚₂=0.1,₁𝑎=𝜇₂=0.1,𝛾=0.25,𝜌=1.25,𝜎=1.05,𝑞=0.35,₁=2.𝜇25,𝑎₂=0.15,𝑏₁=0.75,𝑏₂=1.05,𝑚₁=0.2,𝑚₂=0.1,₁=𝜇₂=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𝜕𝑃𝜕𝑡=𝜕²𝑃𝜕𝑥²𝜕𝑄+𝐹₁(𝑃,𝑄),𝜕𝑡=𝜌𝜕²𝑄𝜕𝑥²+𝜕²𝑅𝜕𝑥²𝜕𝑅+𝐹₂(𝑃,𝑄,𝑅),𝜕𝑡=𝜎𝜕²𝑅𝜕𝑥²+𝛿𝜕²𝑄𝜕𝑥²+𝐹₃(𝑄,𝑅),(4.1) where𝐹₁(𝑃,𝑄)=𝑃(1−𝑃)−𝑎¹𝑃𝑄𝐹𝑚₁+𝑃+𝛾𝑄₂,(𝑃,𝑄,𝑅)=𝑏¹𝑃𝑄𝑚₁+𝑃+𝛾𝑄−𝑎²𝑄²𝑅𝑚²²+𝑄²𝐹−𝜇₁𝑄,₃(𝑄,𝑅)=𝑏²𝑄²𝑅𝑚²²+𝑄²−𝜇₂+𝑞𝑅.(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𝜕𝑃′𝜕𝑡=𝜕²𝑃′𝜕𝑥²𝜕𝑄′+𝑎₁₁𝑃^+𝑎₁₂𝑄^+𝑎₁₃𝑅^,𝜕𝑡=𝜌𝜕²𝑄′𝜕𝑥²+𝜕²𝑅′𝜕𝑥²𝜕𝑅′+𝑎₂₁𝑃^+𝑎₂₂𝑄^+𝑎₂₃𝑅^,𝜕𝑡=𝜎𝜕²𝑅′𝜕𝑥²+𝛿𝜕²𝑄′𝜕𝑥²+𝑎₃₁𝑃^+𝑎₃₂𝑄^+𝑎₃₃𝑅^,(4.5) where𝑎_𝑖𝑗=⎛⎜⎜⎜⎜⎜⎜⎜⎝𝜕𝐹¹𝜕𝜕𝐹𝑃₂′𝜕𝜕𝑄𝐹¹′𝜕𝜕𝑅𝐹¹′𝜕𝜕𝐹𝑃₃′𝜕𝜕𝑄𝐹²′𝜕𝜕𝑅𝐹²′𝜕𝑃′𝜕𝜕𝑄𝐹³′𝜕𝜕𝑅𝐹³′⎞⎟⎟⎟⎟⎟⎟⎟⎠(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|||||||||𝜆+𝑘₂−𝑎−₂₁𝑎₁₁−𝑎₁₂0𝜆+(1+𝜌)𝑘²−𝑎₂₂−𝑎₂₃0−𝑎₃₂𝜆+(𝜌+𝛿)𝑘²−𝑎₃₃|||||||||=0,(4.7) and one obtains the characteristic polynomial𝜆³+𝑑₁₁+𝑑₂₂++𝑑₃₃𝑑𝜆²₁₁+𝑑₁₁𝑑₂₂+𝑑₁₁𝑑₃₃+𝑑₂₂𝑑₃₃−𝑎₁₂𝑎₂₁−𝑎₂₃𝑎₃₂𝜆𝑑₂₂𝑑₃₃−𝑎₂₃𝑎₃₂𝑑₁₁−𝑎₁₂𝑎₂₁𝑑₃₃=0,(4.8) where 𝑑₁₁=𝑘²−𝑎₁₁,𝑑₂₂=(1+𝜌)𝑘²−𝑎₂₂ and 𝑑₃₃=(𝜌+𝛿)𝑘²−𝑎₃₃.

For stability, the following conditions should be satisfied: 𝑑𝑑₁₁+𝑑₂₂+𝑑₃₃>0,(4.9)₁₁𝑑𝑑₂₂𝑑₃₃−𝑎₂₃𝑎₃₂𝑑₁₁−𝑎₁₂𝑎₂₁𝑑₃₃>0,(4.10)₁₁+𝑑₂₂+𝑑>₃₃𝑑𝑑₁₁₁₁𝑑₂₂+𝑑₁₁𝑑₃₃+𝑑₂₂𝑑₃₃−𝑎₁₂𝑎₂₁−𝑎₂₃𝑎₃₂𝑑₂₂𝑑₃₃−𝑎₂₃𝑎₃₂𝑑₁₁−𝑎₁₂𝑎₂₁𝑑₃₃.(4.11)

If diffusion is neglected, that is, 𝑘=0,𝑑₁₁+𝑑₂₂+𝑑₃₃>0⟹𝑎₁₁+𝑎₂₂+𝑎₃₃<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𝐻𝐻𝑘²=𝐴₁𝑘⁶−𝐴₂𝑘⁴+𝐴₃𝑘²+𝐴₄<0,(4.13)_𝑘²=𝐵₁𝑘⁶−𝐵₂𝑘⁴+𝐵₃𝑘²+𝐵₄<0,(4.14)

where 𝐴₁=𝜌+𝜌²,𝐴₂=𝑎₁₁𝜌²+(𝑎₁₁+𝑎₂₂+𝑎₃₃)𝜌+𝑎₃₃,𝐴₃=(𝑎₁₁𝑎₂₂+𝑎₁₁𝑎₃₃−𝑎₁₂𝑎₂₁)𝜌+𝑎₁₁(𝑎₂₂+𝑎₃₃) and  𝐴₄=𝑎₁₁𝑎₂₃𝑎₃₂+𝑎₁₂𝑎₂₁𝑎₃₃−𝑎₁₁𝑎₂₂𝑎₃₃;  while  𝐵₁=(2𝜌+1)(𝜌+2)(𝜌+1),𝐵₂=(2𝜌+1)(2𝜌+3)𝑎₁₁+3(𝜌+1)𝑎²²²+(𝜌+2)(3𝜌+2)𝑎₃₃,  𝐵₃=(2𝜌+1)(𝑎²¹¹−𝑎₂₃𝑎₃₂)+4(𝜌+1)(𝑎₁₁𝑎₂₂+𝑎₁₁𝑎₃₃+𝑎₂₂𝑎₃₃)+(𝜌+2)(𝑎²³³−𝑎₁₂𝑎₂₁)+(𝜌+1)𝑎²²², and  𝐵₄=(𝑎₁₁+𝑎₂₂)(𝑎₁₂𝑎₂₁−𝑎²³³)+(𝑎₂₂+𝑎₃₃)(𝑎₂₃𝑎₃₂−𝑎²¹¹)−(𝑎₁₁+𝑎₃₃)𝑎²²²−2𝑎₁₁𝑎₂₂𝑎₃₃.

The minimum value of 𝐻(𝑘²) occurs at 𝑘²=𝑘^2min=(𝐴₂+(𝐴²²−3𝐴₁𝐴₃))/(3𝐴₁).

Hence for instability to set in, 𝐻𝑘^2min=𝐴²+𝐴²²−3𝐴₁𝐴₃³27𝐴²¹+𝐴³𝐴²+𝐴²²−3𝐴₁𝐴₃3𝐴₁+𝐴₄<𝐴²𝐴²+𝐴²²−3𝐴₁𝐴₃²9𝐴²¹.(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.