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

Authors: J. N. Ndam; J. P. Chollom; T. G. Kassem
Journal: ISRN Applied Mathematics
 (2012)

---

## 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.


---

## Body

## 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)f(P)=cP2a2+P2,
where a is the half-saturation constant and c is the maximal consumption rate of the predator and P is the prey density, while the BD functional response is given by(1.2)f(P,Q)=αPλ+P+γQ,
where Q 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 producerP, called the prey, an intermediate predator Q, and a top predator R, a term adopted from [2]. We assume that the prey P has a logistic growth in the absence of the intermediate predator. It is also assumed that the intermediate predator Q has a BD-functional response, while the top predator R, has a Holling’s type III functional response and that all the three are mobile. Hence we obtain the following model equations:(2.1)∂P∂t=D1∂2P∂x2+rP(1-PK)-a1PQλ+P+γQ,(2.2)∂Q∂t=D2∂2Q∂x2+b1PQλ+P+γQ-a2Q2Rβ2+Q2-μ1Q,(2.3)∂R∂t=D3∂2R∂x2+b2Q2Rβ2+Q2-(μ2+q)R,(2.4)P(x,0)>0,Q(x,0)>0,R(x,0)>0
subject to the reflecting boundary conditions(2.5)D1∂P(0,t)∂x=D1∂P(L,t)∂x=0,D2∂Q(0,t)∂x=D2∂Q(L,t)∂x=0,D3∂R(0,t)∂x=D3∂R(L,t)∂x=0
or the absorbing boundary conditions(2.6)P(0,t)=P(L,t)=0,Q(0,t)=Q(L,t)=0,R(0,t)=R(L,t)=0,
where r>0 is the intrinsic growth rate of the prey and K>0 is its carrying capacity and D1, D2, and D3 are the constant diffusion coefficients of the prey, the intermediate predator, and the top-predator, respectively; μ1>0, μ2>0, and q>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 ai, bj, β, λ, 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:(2.7)t=t*r,P=KP*,Q=KQ*,R=KR*,x=Lx*,a1=ra1*,a2=ra2*,b1=rb1*,b2=rb2*,μ1=rμ1*,μ2=rμ2*,q=rq*,
where the asteriated variables are nondimensional variables. We obtain after dropping asterisks the scaled system(2.8)∂P∂t=∂2P∂x2+P(1-P)-a1PQm1+P+γQ,(2.9)∂Q∂t=ρ∂2Q∂x2+b1PQm1+P+γQ-a2Q2Rm22+Q2-μ1Q,(2.10)∂R∂t=σ∂2R∂x2+b2Q2Rm22+Q2-(μ2+q)R,(2.11)P(x,0)>0,Q(x,0)>0,R(x,0)>0,(2.12)∂P(0,t)∂x=∂P(L,t)∂x=0,ρ∂Q(0,t)∂x=ρ∂Q(L,t)∂x=0,σ∂R(0,t)∂x=σ∂R(L,t)∂x=0,
where we have taken L=D1/r, ρ=D3/D1,σ=D2/D1, m1=λ/K, and m2=β/K.


### 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 calledcross-diffusion. We introduce in this model cross-diffusion terms with constant coefficients D12 and D21. The scaled model equations then take the form(2.13)∂P∂t=∂2P∂x2+P(1-P)-a1PQm1+P+γQ,(2.14)∂Q∂t=ρ∂2Q∂x2+∂2R∂x2+b1PQm1+P+γQ-a2Q2Rm22+Q2-μ1Q,(2.15)∂R∂t=σ∂2R∂x2+δ∂2Q∂x2+b2Q2Rm22+Q2-(μ2+q)R,
where δ=D21/D12, subject to the boundary conditions (2.6) and (2.12).


## 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 calledcross-diffusion. We introduce in this model cross-diffusion terms with constant coefficients D12 and D21. The scaled model equations then take the form(2.13)∂P∂t=∂2P∂x2+P(1-P)-a1PQm1+P+γQ,(2.14)∂Q∂t=ρ∂2Q∂x2+∂2R∂x2+b1PQm1+P+γQ-a2Q2Rm22+Q2-μ1Q,(2.15)∂R∂t=σ∂2R∂x2+δ∂2Q∂x2+b2Q2Rm22+Q2-(μ2+q)R,
where δ=D21/D12, 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(3.1)a1=1.2,a2=0.25,b1=0.5,b2=0.75,m1=0.2,m2=0.1,μ1=μ2=0.1,γ=0.25,ρ=1.25,σ=1.05,q=0.35,a1=2.25,a2=0.15,b1=0.75,b2=1.05,m1=0.2,m2=0.1,μ1=μ2=0.1,γ=0.5,ρ=1.10,σ=1.5,q=0.5,
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.Figures6 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)∂P∂t=∂2P∂x2+F1(P,Q),∂Q∂t=ρ∂2Q∂x2+∂2R∂x2+F2(P,Q,R),∂R∂t=σ∂2R∂x2+δ∂2Q∂x2+F3(Q,R),
where(4.2)F1(P,Q)=P(1-P)-a1PQm1+P+γQ,F2(P,Q,R)=b1PQm1+P+γQ-a2Q2Rm22+Q2-μ1Q,F3(Q,R)=b2Q2Rm22+Q2-(μ2+q)R.
Suppose that there exists a spatially uniform steady state (Pc,Qc,Rc), such that(4.3)Fj(Pc,Qc,Rc)=0,j=1,2,3.
To determine the stability of the uniform steady state, we perturb the population densities thus(4.4)P(x,t)=Pc+P′(x,t),Q(x,t)=Qc+Q′(x,t),R(x,t)=Rc+R′(x,t).
Hence the linearised system becomes(4.5)∂P′∂t=∂2P′∂x2+a11P′+a12Q′+a13R′,∂Q′∂t=ρ∂2Q′∂x2+∂2R′∂x2+a21P′+a22Q′+a23R′,∂R′∂t=σ∂2R′∂x2+δ∂2Q′∂x2+a31P′+a32Q′+a33R′,
where(4.6)(aij)=(∂F1∂P′∂F1∂Q′∂F1∂R′∂F2∂P′∂F2∂Q′∂F2∂R′∂F3∂P′∂F3∂Q′∂F3∂R′)
is the Jacobian matrix for the system evaluated at the uniform equilibrium point (Pc,Qc,Rc). In order to examine the linear stability of the system (4.5), we use the normal modes approach by assuming solutions of the form P′~e(λt+ikx),Q′~e(λt+ikx), and R′~e(λt+ikx), where λ and k are constants. Thus the eigenvalues satisfy the equation(4.7)|λ+k2-a11-a120-a21λ+(1+ρ)k2-a22-a230-a32λ+(ρ+δ)k2-a33|=0,
and one obtains the characteristic polynomial(4.8)λ3+(d11+d22+d33)λ2+(d11d22+d11d33+d22d33-a12a21-a23a32)λ+(d11d22d33-a23a32d11-a12a21d33)=0,
where d11=k2-a11,d22=(1+ρ)k2-a22 and d33=(ρ+δ)k2-a33.For stability, the following conditions should be satisfied:(4.9)d11+d22+d33>0,(4.10)d11d22d33-a23a32d11-a12a21d33>0,(4.11)(d11+d22+d33)(d11d22+d11d33+d22d33-a12a21-a23a32)>d11d22d33-a23a32d11-a12a21d33.If diffusion is neglected, that is,k=0,(4.12)d11+d22+d33>0⟹a11+a22+a33<0.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(4.13)H(k2)=A1k6-A2k4+A3k2+A4<0,(4.14)H′(k2)=B1k6-B2k4+B3k2+B4<0,whereA1=ρ+ρ2,A2=a11ρ2+(a11+a22+a33)ρ+a33,A3=(a11a22+a11a33-a12a21)ρ+a11(a22+a33) and  A4=a11a23a32+a12a21a33-a11a22a33;  while  B1=(2ρ+1)(ρ+2)(ρ+1),B2=(2ρ+1)(2ρ+3)a11+3(ρ+1)a222+(ρ+2)(3ρ+2)a33,  B3=(2ρ+1)(a112-a23a32)+4(ρ+1)(a11a22+a11a33+a22a33)+(ρ+2)(a332-a12a21)+(ρ+1)a222, and  B4=(a11+a22)(a12a21-a332)+(a22+a33)(a23a32-a112)-(a11+a33)a222-2a11a22a33.The minimum value ofH(k2) occurs at k2=kmin2=(A2+(A22-3A1A3))/(3A1).Hence for instability to set in,(4.15)H(kmin⁡2)=(A2+A22-3A1A3)327A12+A3(A2+A22-3A1A3)3A1+A4<A2(A2+A22-3A1A3)29A12.Similar result is obtained forH′(kmin2) in terms of the Bi,i=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 (Figure1). 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 Figures2 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.



---

*Source: 391547-2012-03-08.xml*