Prey, predator and super-predator model with disease in the super-predator

ARTICLE IN PRESS

JID: AMC

[m3Gsc;October 28, 2016;20:34]

Applied Mathematics and Computation 0 0 0 (2016) 1–23

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Prey, predator and super-predator model with disease in the super-predator W. Mbava a,∗, J.Y.T. Mugisha a,b, J.W. Gonsalves a a

Department of Mathematics & Applied Mathematics, Nelson Mandela Metropolitan University, P.O. Box 77000, Port Elizabeth 6031, South Africa b Department of Mathematics, Makerere University, P.O. Box 7062, Kampala, Uganda

a r t i c l e

i n f o

Keywords: Predator–prey system Eco-epidemiology Stability Infectious disease

a b s t r a c t The dynamics of a predator–prey model with disease in super-predator are investigated. The predator is under immense competition from the super-predator and is also facing extinction. The disease is considered as biological control to allow the predator population to regain from a low number. The results highlight that in the absence of additional mortality on predator by super-predator, the predator population survives extinction. At current levels of disease incidence, the super-predator population is wiped out by the disease. However, the super-predator population survives extinction if the disease incidence rate is low. Persistence of all populations is possible in the case of low disease incidence rate and no additional mortality imparted on predator. Furthermore, a two-species subsystem, prey and predator, is considered as a special case to determine the effect of super-predator removal from the system, on the survival of the predator. This is treated as a contrasting case of the smaller parks. The results show that the predator population thrives well in the total absence of its main competitor, with its population rising to at least twice the initial value. © 2016 Elsevier Inc. All rights reserved.

1. Introduction The sizes of species populations are important in ecological studies. They are influenced by ecological and epidemiological factors. The ecological factors include species interactions in the form of competition and predation. The epidemiological factors include the spread of infectious diseases [1]. The study of transmissible diseases within an ecological setting is gaining momentum [2]. It is becoming biologically relevant to include the effects of diseases in studies on the behaviour of dynamical ecological systems [3]. The effect of disease in prey on prey–predator systems has been studied by several researchers. For most such models, the paramount assumption is that predation favours infected rather than sound prey [2]. Lu et al. [4] proposed a model for hunting strategies of predators in susceptible and infected prey. Their results show that predation is a function of the strength of the disease. When the disease is slightly endemic, predators target the susceptible prey, but switch to infected prey when the disease becomes heavily endemic. Chakraborty et al. [5] studied a prey–predator eco-epidemiological system on disease persistence and extension perspective. They outlined a criterion for describing uniformly strong persistence of the system. The disease affected the prey. Their results indicate that the disease incidence rate has a potential to ∗ Corresponding author at: Department of Mathematics & Applied Mathematics, Nelson Mandela Metropolitan University, P.O. Box 770 0 0, Port Elizabeth 6031, South Africa. E-mail addresses: [email protected], [email protected] (W. Mbava).

http://dx.doi.org/10.1016/j.amc.2016.10.034 0 096-30 03/© 2016 Elsevier Inc. All rights reserved.

Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

JID: AMC 2

ARTICLE IN PRESS

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

destabilise a once stable system. Sahoo and Poria [6] studied the diseased prey–predator model with general Holling type interactions. Their model involved prey, an intermediate predator, and a top predator. They investigated different types of general Holling interactions. Their findings indicate that by properly choosing interaction functions, a diseased system can be transformed into a disease-free system, thereby controlling the disease. Hethcote et al. [7] studied a predator–prey model with Susceptible–Infective–Susceptible (SIS) disease dynamics in the prey. The infected prey was more vulnerable to predation. They showed that some parameter values increased the risk of predation of infected prey, thereby ensuring that the predator population persists. Otherwise the predators would become extinct in the absence of disease. Furthermore, they showed that the increased predation of vulnerable infected prey would allow the disease to die out, apart from the disease remaining endemic in the absence of predators. Mukhopadhyay and Bhattacharyya [3] studied the dynamics of a delay-diffusion prey–predator model with disease in prey. For their basic model, they established that persistence of the disease depended on the predator death rate and the basic reproduction number. Upon including delay attributed to gestation of the predator, they found out that the predator death rate, basic reproduction number, and equilibrium density of susceptible prey together shaped the dynamical behaviour of the system. On the role of diffusion in the delayed model, they deduced that diffusivity coefficients for susceptible and infected prey, together with the perturbation wave number of the general solution, determined the dynamical behaviour of the system. Haque et al. [8] investigated an eco-epidemiological predator–prey model with disease in prey. They showed that a virulent disease in the prey may allow predators to escape extinction but destabilises a once stable system. Xiao and Chen [9] proposed a predator–prey model with disease in prey. Their model shows that the introduction of a time delay in the coefficient of converting prey into predators has both stabilising and destabilising effects on the positive steady state. Sahoo and Poria [10] studied the effects of additional food on an eco-epidemic model with time delay on infection. The prey had a parasitic infection. Additional food was provided to predators as a means of controlling its population, and ultimately the parasitic infection prey. Their results show that the system becomes disease free under provision of additional food to predators. Predator population increases thereby enhancing the consumption of diseased prey. Sahoo [11] investigated the role of additional food provided to predators in an eco-epidemiological system with disease in prey. An optimal control problem was formulated and solved to determine the control of disease. The results show existence of a critical infection rate above which disease free system is not reachable without additional food. However, additional food enables a disease free system to be attained up to a certain infection level. Other researchers consider the situation where the disease spreads among the predator population. The predator–prey interactions are extended to include disease in the predator species. Haque [2] studied the predator–prey model with SIS parasitic infection spreading through the predator species only. It was shown that infection in the predator species may save the prey from extinction even if the basic reproduction number for the prey to be able to invade the predator-only equilibrium, was less than one. Pal et al. [12] investigated a predator–prey model with disease present in predator species only. They showed that for some values of the predation rate all species could be saved from extinction, and that the disease did not propagate in the predator population. Han et al. [13] studied four predator prey models in which disease spreads in both the prey and predator. The disease transmission involved both mass action and standard incidence. They showed that when the disease persists in the prey population and the predators feed sufficiently to survive, the disease also persists in the predator population. In this study, the model proposed is an extension of Haque [2] and Venturino [14] models to include an infectious disease spreading through a super-predator species. The study area is Kruger National Park (KNP), South Africa’s largest wildlife reserve with an area of approximately 20,0 0 0 km2 . The park consists largely of woodland savanna. The prey, predator and super-predator species are the impala, cheetah and lion, respectively. It is considered that the cheetah feeds only upon the impala, the lion feeds upon the impala but kills the cheetah to reduce competition. The impala is the most preferred prey species amongst a host of species consumed by the cheetah [15,16]. However, it is the third dominant prey species consumed by the lion after zebra and wildebeest in the KNP [17]. Despite the impala occurring in high numbers, the cheetah population remains very low. Moreover, the IUCN (World Conservation Union) Red Data Books listed the cheetah species as vulnerable in South Africa [18]. Bovine Tuberculosis (BTB) disease, caused by Mycobacterium bovis, is present in the lion population of KNP [19,20]. The initial reported outbreak of BTB in lions happened in 1995, and was believed to have originated from feeding on buffalo infected by the same disease [19,21]. An 80% disease prevalence in the lion population occupying the southern part of the KNP was reported in the year 20 0 0 [19,22]. Since lions are the apex predators, it has been hypothesized that there may be a change in the lion population, which will ultimately result in changes in the dynamics of other predators, based on their interactions [23]. The aim of this study is to determine whether or not the presence of the disease in the lion population acts as biological control in weakening the lion and indirectly improving the population of the cheetah. With the cheetah confronted with major competition from the lion, as well as extinction, the study also determines the extent to which the presence of the lion affects the cheetah population. 2. Model formulation 2.1. Assumptions The model consists of three populations, the impala, cheetah and lion species whose populations are denoted by U(t), V(t) and N(t), respectively. The following assumptions are considered in formulating the model. Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

[m3Gsc;October 28, 2016;20:34] 3

Table 1 Parameters table for systems (2.1) and (2.2). Parameter

Description

r K a b c e f l m n p

Intrinsic growth rate of the prey Environmental carrying capacity of the prey Capture rate of the prey by the predator Capture rate of the prey by the super-predator Efficiency of infected super-predators to capture prey Handling time of prey by predator Handling time of prey by super-predator Conversion rate of prey biomass into new predators Conversion rate of prey biomass into new susceptible super-predators Conversion rate of prey biomass into new infected super-predators Mortality rate of predator by super-predator Disease standard incidence Natural mortality rate of predators Natural mortality rate of super-predators Disease-induced mortality rate of infected super-predators

β μ

ν δ

(A1) In the presence of disease, the lion population is divided into two classes, namely susceptible lion and infected lion denoted by W(t) and B(t) respectively. Therefore at time t the total lion population is N (t ) = W (t ) + B(t ). (A2) The disease spreads among the lion population only through in-group and out-group interactions and is not genetically inherited. The infected population do not recover or become immune. The incidence function β (N) is assumed to be the nonlinear function β (N ) ≡ N1 , the standard incidence. (A3) The cheetah and lion consume the impala, although with possibly different search efficiencies, denoted by a and b. (A4) The cheetah and lion have natural death rates μ and ν , respectively. Furthermore, infected lion have disease-induced death rate δ . The susceptible lion kills off the cheetah at a rate p. Infected lions are considered too weak to kill cheetah.

2.2. The model In the super predator–predator–prey model, the cheetah and lion consume the impala according to the nonlinear Holling type II functional response. The Holling Type II response is the most common type of functional response for predator species [24]. The impala are easier to catch. A solitary lion or cheetah can catch an impala as the impala do not exhibit group defence. On the basis of the assumptions above the following model is proposed





U aUV bUN dU = rU 1 − − − , dt K 1 + eU 1 + fU dV lUV = − pV N − μV, dt 1 + eU dN mUN = − ν N, dt 1 + fU

(2.1)

where U, V, and N represent prey, predator and super-predator populations. The Susceptible–Infective (SI) dynamics are incorporated into the model (2.1) for the disease in super-predator using the model proposed by Venturino [14]. The model (2.1) becomes

dU dt dV dt dW dt dB dt dN dt



= rU 1 − lUV 1 + eU mUW = 1 + fU nUB = 1 + fU mUW = 1 + fU =

U K



−

aUV bUW cbUB − − , 1 + eU 1 + fU 1 + fU

− pV W − μV,

(2.2)

βW B

− νW, N βW B + − ν B − δ B, N nUB + − δ B − ν N, 1 + fU −

where n = cm, W and B are the susceptible and infected super-predator respectively, and N = W + B. The parameters used in the model are explained in Table 1. Let S and I denote the fractions of susceptible and infected super-predators defined Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC 4

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

by S =

W B and I = such that S + I = 1. Since N N

I =

 B  N

=





B 1  B − N , N N

S =

 W  N

=





1  W  B − N . N N

S and I can be expressed in terms of W and B , respectively. After some calculations and dropping the redundant equation involving the susceptible super-predator S , the system (2.2) reduces to

dU dt dV dt dN dt dI dt



= U r 1− =V

 lU

U K

1 + eU



−



aV bN[(1 − I ) + cI] − = F1 (U, V, N, I ), 1 + eU 1 + fU



− p(1 − I )N − μ = F2 (U, V, N, I ),

 [m(1 − I ) + nI]U

(2.3)



− δ I − ν = F3 (U, V, N, I ), 1 + fU   (m − n )U = I (1 − I ) β − − δ = F4 (U, V, N, I ), =N

1 + fU

with the initial conditions 0 ≤ U0 ≤ K, 0 ≤ V0 ≤ A, 0 ≤ N0 ≤ N1 and 0 ≤ I0 ≤ 1. Thus, the system (2.3) is studied in the region  = {(U, V, N, I ) ∈ R4+ : 0 ≤ U ≤ K, 0 ≤ V ≤ A, 0 ≤ N ≤ N1 , 0 ≤ I ≤ 1}. 3. Model analysis 3.1. Boundedness All the parameters of system (2.3) are non-negative, and so the corresponding right-hand side of the system is a smooth function of the variables (U, V, N, I) in the region . It follows that local existence and uniqueness properties hold for the solution of the system. Proposition 1. R4+ is an invariant set. Proof. The system (2.3) is homogeneous. It follows that the coordinate planes U = 0, V = 0, N = 0 and I = 0 are solutions for it. By the existence and uniqueness theorem, any trajectory emanating from R4+ stays there and cannot cross the coordinate planes. Thus, R4+ is an invariant set.  Proposition 2.  = {(U, V, N, I ) ∈ R4+ : 0 ≤ U ≤ K, 0 ≤ V ≤ A, 0 ≤ N ≤ N1 , 0 ≤ I ≤ 1} in invariant under the flow (2.3).

) = U Kr (K − U ) with solution U (t ) ≤ D+KD e−rt lU and lim sup U (t ) ≤ K by the comparison principle for ODEs. From the second equation, dV ≤ V ( − μ ) ≤ V ( l − μ ) if e 1+eU dt Proof. Consider the system (2.3). The first equation implies that

dU dt

≤ Ur (1 −

U K

t→∞

> 1. The corresponding solution is V (t ) ≤ Ae(l−μ )t and lim sup V (t ) ≤ A if l < μ. The third equation implies that −I )+nI]U N[ [m(11+ fU

− ν ] ≤ N (mK + nK − ν ) with solution N (t ) ≤

lows from the fourth equation that

dI dt

t→∞ N1 e[K (m+n )−ν ]t

≤ β I (1 − I ) with solution I (t ) ≤

dN dt

≤

and lim sup N (t ) ≤ N1 when (m + n )K < ν . It fol-

D D + e −β t

t→∞

and lim sup I (t ) ≤ 1. t→∞

Any solution with initial values U0 , V0 , N0 and I0 eventually approaches this region.



¯ are uniformly bounded if l < a, n < m, m < b and c < 1. Proposition 3. All solutions of system (2.3) which initiate in  Proof. Define the function

 = U + V + N.

(3.1)

Differentiating  with respect to time along the solutions of (2.3) gives





U aUV bNU[(1 − I ) + cI] d lUV [m(1 − I ) + nI]UN = Ur 1 − − − + − p(1 − I )NV − μV + − δ IN − ν N] dt K 1 + eU 1 + fU 1 + eU 1 + fU   U ≤ Ur 1 − − μV − ν N K   rU = U r+φ− − (μ − φ )V − (ν − φ )N − φ . K Thus,

d r K (r + φ ) + φ  = − [U 2 − U] − (μ − φ )V − (ν − φ )N dt K r   r K ( r + φ ) 2 K ( r + φ )2 =− U− + − (μ − φ )V − (ν − φ )N K 2r 4r Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

≤

5

K ( r + φ )2 − (μ − φ )V − (ν − φ )N. 4r

Choosing φ to be such that φ ≤ min {μ, ν } the right-hand side will be bounded. Let ϕ be such that

K ( r + φ )2 d + φ ≤ = ϕ. dt 4r

(3.2)

A solution of this is

(t ) ≤ Ce−φt +

ϕ φ

and

(t ) ≤ (0 )e−φt +

 ϕ ϕ (1 − e−φt ) ≤ max (0 ), . φ φ

ϕ Moreover, lim sup (t ) ≤ φ as t → ∞ independent of initial conditions.



3.2. Steady states Since the focus is on the growth of animal species, there is need for the steady states of the system to satisfy conditions for non-negativity. Furthermore, it is realised that the predators cannot survive in the complete absence of their prey. As such, the steady states E(0, V, N, I), E(0, 0, N, I), E(0, 0, 0, I), E(0, V, 0, 0), E(0, 0, N, 0) and E(0, 0, 0, I) are not feasible. However, the steady state E(0, 0, 0, 0) always exist. The steady states in which prey is present are studied only. The system (2.3) has four biologically feasible nonnegative steady states. (i) E1 (U1 , 0, N1 , 0) where U1 =

ν

m− f ν

and N1 =

rm (K m−ν − f K ν ) . bK ( f ν −m )2

Prey and super-predators coexist in the absence of disease

and predators. The equilibrium point E1 exists if m > fν , growth function of super-predators exceed their death, and K (m − f ν ) > ν or K > m−ν f ν , prey grow to a population level which is less than their carrying capacity.

U bN2 −mμ−eμν , N2 = l ν +pf(μν and V2 = Da , where D = (1 + eU2 )[r (1 − K2 ) − 1+ fU . Prey, m + eν − f ν ) 2 predator and super-predator coexist in the absence of disease. The equilibrium point E2 exists if m > fν and ν (l − U2 μ eμ ) > μ(m − f ν ), that is, m−ν f ν > l−e μ , the growth ratio of super-predators exceed that of predators, r (1 − K ) >

(ii) E2 (U2 , V2 , N2 , 0) where U2 =

ν

m− f ν

the growth function of prey exceeds their predation rate by super-predators, and m + eν − f ν = m − f ν + eν > 0. −δ (iii) E3 (U3 , 0, N3 , I3 ), where U3 = m−nβ−−f βδ + f δ = (m−n β )− f (β −δ ) , bN2 , 1+ fU2

I3 =

mβ −mδ −mν +nν β (m−n )

C=

=

m (β −δ )−ν (m−n ) , β (m−n )



and N3 =

C G,

where



r (m − n ) U3 − 1 [U3 [(m − n ) + f δ ] + δ ] < 0 b K

as long as U3 < K, and

G = (n − m )(ν + cβ I ) + (nβ − mδ ) < 0 if nβ − mδ < 0 or βδ < m n , that is, the ratio of disease incidence to disease induced death is less than that of the conversion rate of prey biomass into new susceptible super-predators to infected ones. Hence, N3 > 0. Prey and superpredator coexist in the presence of disease but in the absence of the predator. The equilibrium point E3 exists if m > n, conversion rate of prey biomass into new super-predators for susceptible super-predators exceeds that of infected ones, β > δ , disease standard incidence exceeds disease-induced mortality rate of infected super-predators, and (m − n ) > f (β − δ ), m(β − δ ) > ν (m − n ), that is, m ν (β − δ ) > (m − n ) > f (β − δ ) or m − f ν > 0 (condition for existence of E2 ) and N3 > 0. lU

4 ) β (m−n )(μ− 1+eU

−δ m (β −δ )−ν (m−n ) (iv) E4 (U4 , V4 , N4 , I4 ) given by U4 = (m−n β , N4 = p(nβ −mδ +(n−m )4ν ) and V4 = 1a (r (1 + eU4 )(1 − )− f (β −δ ) , I4 = β (m−n ) U4 bN4 K ) + beβ N4U4 I4 + 1+ fU (I4 − β (1 + cI4 )(1 − eU4 )). U4 > 0 and I4 > 0 (conditions for existence of E3 ). Prey, predator 4

and super-predator coexist in the presence of disease. This equilibrium point contrasts E2 and exposes the effects of lU4 the presence of disease on the system dynamics. The equilibrium point E4 exists if μ < 1+eU , predator’s death rate 4 is less than its growth function, and (nβ − mδ ) + (n − m )ν < 0 from condition of existence of E3 , hence, N4 > 0; and U4 < K and I4 > β (1 + cI4 )(1 − eU4 ), in which case V4 > 0. Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC 6

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

3.3. Analysis of steady states The Jacobian matrix for the system (2.3) is given by

⎡ ⎢ ⎢

lV

lU 1+eU

J=⎢ ⎢ (m−(m−n)I )N (1+eU )2

⎣

1+cI−I ) − bU (1+ fU

−aU 1+eU

j11

− N p( 1 − I ) − μ

(1+ fU )2 m−n )I − ((1+ fU )2

−pV (1 − I )

0

U (m−mI+nI ) 1+ fU

0

0

⎤

bUN 1+ fU

⎥ ⎥ ⎥, (m−n )UN − 1+ fU − δ N⎥ ⎦ β − (m1+−nfU)U − δ pNV

− δI − ν

(1+cI−I ) − (1+aV − bN(1+ . At any steady state solution, the Jacobian matrix is computed. Let Jk = J denote eU )2 fU )2 [k] the Jacobian evaluated at Ek and ji j = bi j , i = 1, 2, 3, 4, j = 1, 2, 3, 4, k = 1, 2, 3, 4, the corresponding entries. 2rU K

where j11 = r −

3.3.1. Local and global stability of the steady state E1 (K m−ν − f K ν ) For the equilibrium point E1 ( m−ν f ν , 0, rmbK , 0 ), the Jacobian matrix is given by ( f ν −m )2

⎡

r−

2rU1 K

− (1+bNfU1 )2 1 0

⎢ ⎢ J1 = ⎢ ⎢ − f mU1 N1 + ⎣ 1+ fU1

− 1+aUeU1 1 −μ − pN1 +

mN1 1+ fU1

0

1 − 1+bUfU

bU1 N1 1+ fU1

⎥ ⎥ ⎥. (m−n )U1 N1 ⎥ −δ N1 − 1+ fU ⎦ 1 β − δ − (m1+−nfU)U1 1

1

lU1 1+eU1

0

0

mU1 1+ fU1

0

0

⎤

0

−ν

The eigenvalues of J1 are





−b1 + b21 −4b0 b2 −b1 − b21 −4b0 b2 [1] [1] [1] ν− f Kν ) λ[1] = β − δ − ν + nmν , λ2 = −μ − prmbK(K(mf ν−−m + m+(leν− f )ν ; λ3 = and λ4 = , for which 1 2b0 2b0 )2 [1] [1] 2 λ3 and λ4 are the roots of the equation b0 λ + b1 λ + b2 = 0 where b0 = Km(m − f ν ), b1 = rν [(m + f ν ) − f K (m − f ν )] and b2 = rν [K (m + f ν )2 + f 2 ν 2 ) − ν (m − f ν ) + 2erN13 (l − eμ ) − aK μV1 . Note that b0 > 0, that is, condition for existence of 2 +fν equilibrium point E1 . Thus, b1 > 0 if m > f K. b2 > 0 if ((mm+−ffνν)) > Kν and 2erN13 (l − eμ ) > aK μV3 . m− f ν

Lemma 1. The system (2.3) is locally asymptotically stable about E1 if (a)

ν− f Kν ) lν < μ + prmbK(K(mf ν−−m , that is, the sum m+ ( e− f )ν )2 nν β + m < δ + ν, that is, mβ + nν < m(δ + ν ),

(b) (c) b1 > 0, b2 > 0, m > fν and K >

ν

m− f ν

of the predator’s killing and death rates exceeds their growth rate,

.

Now, the global stability of E1 is discussed whenever it exists, using the Bendixson–Dulac criterion, Theorem 5, [25] and the Poincar e´ –Bendixson Theorem, Theorem 6, [26]. Proposition 4. E1 is globally asymptotically stable if 2ν > mK. Proof. Consider subsystem (ii) of system (2.3)





U bN dU = Ur 1 − − , dt K 1 + fU dN mNU = − ν N. dt 1 + fU Let (U, N ) = It follows that

1 UN ,

h1 (U, N ) = Ur (1 −

(3.3) U K

)−

bUN , 1+ fU

h2 (U, N ) =

mUN 1+ fU

− μN. Now (U, N) > 0 in the interior of the U − N plane.

∂ ( h1 ) ∂ ( h2 ) + ∂U ∂N     ∂  m U b ∂ r ν = 1− − + − ∂U N K 1 + fU ∂ N 1 + fU U b f KN − r (1 + fU )2 = . KN (1 + fU )2

(U, N ) =

Now  is not identically zero in the positive quadrant of the U − N plane. Thus, (U, N) < 0 and does not change sign if

b f KN − r (1 + fU )2 = b f KN − r − 2r fU − f 2U 2 = −r − (2rU − bKN ) − f 2U 2 < 0, 2ν −mK )+ν rm that is, if 2rU − bKN > 0. This means that at E1 , we have 2rU1 − bKN1 = r (m− f ν()( > 0 if and only if 2ν > mK. Thus, m− f ν )2

¯ . Since E1 is locally asymptotically stable, the Poincar e´ –Bendixson Theorem subsystem (3.3) does not have a limit cycle in  ¯.  (Theorem 6) and Proposition 3 imply that the equilibrium point E1 is globally asymptotically stable in the region  Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

7

3.3.2. Local and global stability of the steady state E2 For the equilibrium point E2 (U2 , V2 , N2 , 0) the Jacobian matrix is given by

⎡

r−

⎢ ⎢ J2 = ⎢ ⎢ ⎣

2rU2 K

− (1+bNfU2 )2 − (1+aVeU2 )2 2 2

− 1+aUeU2 2

lV2

lU2 1+eU2

(1+eU2 )2

2 − 1+bUfU

⎥ ⎥ ⎥. (m−n )U2 N2 ⎥ −δ N2 − 1+ fU ⎦ 2 β − δ − (m1+−nfU)U2 2

2

− μ − pN2

−pV2

pN2V2

−ν

mN2 1+ fU2

0

mU2 1+ fU2

0

0

0

⎤

bU2 N2 1+ fU2

The local stability of the disease-free equilibrium at E2 is deduced from the eigenvalues of the Jacobian matrix J2 . From (m−n )U the fourth row this matrix has the eigenvalue β − δ − 1+ fU 2 , and the remaining eigenvalues are derived from the 3 × 3 2

Jacobian matrix given by

⎡

⎤

⎡

⎥

⎢

[2] j11

[2] d11

[2] d12

[2] d13

[2] J2 = ⎣d21

[2] d22

[2] d23 ⎦ = ⎣ (1+eU2 2 )2

⎢

[2] d31

[2]

[2]

[2]

lV

[2]

[2]

[2]

⎥ ⎦,

2

−pV2 mU2 1+ fU2

0

[2]

⎤

2 − 1+bUfU

− μ − pN2

lU2 1+eU2

mN2 1+ fU2

[2] d33

0

− 1+aUeU2 2

[2]

[2]

(3.4)

−ν

[2]

[2]

[2]

[2]

[2]

[2]

[2]

[2]

where d11 = j11 , d12 = j12 < 0, d13 = j13 < 0, d21 = j21 > 0, d22 = j22 , d23 = j23 < 0, d31 = j31 > 0, d32 = 0, d33 = j33 . The characteristic polynomial of the Jacobian of the system evaluated at this point is given by

λ3 + a1 λ2 + a2 λ + a3 = 0,

(3.5)

where [2] [2] [2] a1 = −(d11 + d22 + d33 ) = −Tr(J2 ), [2] [2] [2] [2] [2] [2] [2] [2] [2] [2] a2 = d11 d22 + d11 d33 + d22 d33 − d12 d21 − d13 d31 , [2] [2] [2] [2] [2] [2] [2] [2] [2] [2] [2] [2] a3 = −d12 d23 d31 + d13 d22 d31 + d12 d21 d33 − d11 d22 d33 = −Det J2 .

[2]

[2]

[2]

Now, if d11 < 0, d22 < 0, d33 < 0 then a1 > 0 and a2 > 0. a3 > 0 if [2] [2] [2]

[2] [2] [2]

[2] [2] [2]

[2] [2] [2]

d13 d22 d31 + d12 d21 d33 − d11 d22 d33 > d12 d23 d31 . Furthermore, [2] [2] [2] 2[2] [2] [2] [2] [2] [2] 2[2] [2] [2] [2] [2] [4] [4] a1 a2 − a3 = d11 d12 d21 − d11 d22 + d12 d21 d22 − d11 d22 + d11 d13 d31 + d12 d23 d31 2[2] [2] [2] [2] [2] 2[2] [2] [2] [2] [2] [2] 2[2] [4] 2[4] − d11 d33 − 2d11 d22 d33 − d22 d33 + d13 d31 d33 − d11 d33 − d22 d33

> 0. Application of the Routh–Hurwitz criterion [27] leads to the conclusion that the characteristic equation (3.5) has negative real parts. Thus, the following result can be stated. Lemma 2. The equilibrium point E2 is locally asymptotically stable if β > δ +

(m − n )U2 1 + fU2

[2]

[2]

[2]

, d11 < 0, d22 < 0 and d33 < 0.

To analyse global stability of the equilibrium point E2 , the last equation in system equation (2.3) is dropped to obtain a reduced system since the equilibrium point does not depend on I. The Jacobian matrix of the reduced system is given by (3.4). The reduced system is not competitive and so the geometric approach to global stability problems method is used. The framework for this method is developed in the papers of Smith [28] and Li and Muldowney [29] and was used by Buonomo and Lacitignola [30], Li et al. [31] and Tian and Wang [32]. This framework is used to prove the global stability of the equilibrium point E2 . The second compound matrix of the reduced system is

⎡



⎢

J2[2] = ⎣

[2] j11 +

lU2 1+eU2

− μ − pN2

[2] j11 +

mU2 1+ fU2 lV2 (1+eU2 )2

0 mN2 − 1+ fU

2

The matrix function P is set by

⎤

bU2 1+ fU2

−pV2 −ν

⎥ ⎦.

− 1+aUeU2 2 lU2 1+eU2

− μ − pN2 +

mU2 1+ fU2

−ν



V2 V2 P (U2 , V2 , N2 ) = diag 1, , . N2 N2 Then Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC 8

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

⎡

0

⎢ PF = ⎢0 ⎣

V2 N2



0 V

0 and

−

2

V2

N2 N2

V2 N2

0

⎡

0

0

V

PF P −1 = ⎣0

−

2

V2

0 Now

PF J2[2] P −1 = ⎣ 

[2] j11 +



0 V2 V2

0

2

N2

V2 V2

−

P −1

N2 N2

N2 N2

−

1 = ⎣0 0

⎤

0

0 0⎦

N2 V2

N2 V2

0

⎦.

− μ − pN2

lU2 1+eU2

⎡

⎥ ⎥, ⎦ 

⎤

0

N

0

⎡

⎤

0

 

−pN2 [2] j11 +

mU2 1+ fU2 lV2 (1+eU2 )2

0 V2 2 − (1+mN fU2 )2 N2

−ν lU2 1+eU2

bU2 N2 1+ fU2 V2 − 1+aUeU2 2

⎤ ⎦.

− μ − pN2 +

mU2 1+ fU2

−ν

 [2]

The matrix Q = PF P −1 + PF J2 P −1 can be written in block form:





Q11 Q= Q21 [2]

with Q11 = j11 +

Q12

Q12 , Q22 lU2 − μ − pN2 , 1 + eU2



bU2 N2 1+ fU2 V2

= −pN2

 Q22 =



,

 Q21 =



0 2 − (1+mN fU )2 2

[2] j11 +

mU2 1+ fU2

−ν +

V2 V2

−

N2 N2

lV2

,



− 1+aUeU2 2

lU2 1+eU2

(1+eU2 )2

V2 N2

− μ − pN2 +

mU2 1+ fU2

−ν +

V2 V2

−

N2 N2

.

A norm in R3 is defined as |u, v, n| = max{|u|, |v| + |n|} for any vector (u, v, n ) ∈ R3 . Let m denote the Lozinski ˘i measure with respect to this norm. Then

m(Q ) ≤ sup{g1 , g2 }

(3.6)

with g1 = m1 (Q11 ) + |Q12 | and g2 = |Q21 | + m1 (Q22 ), where |Q12 | and |Q21 | are matrix norms induced by the L1 vector norm, and m1 denotes the Lozinski ˘i measure with respect to the L1 norm. More specifically, lU2 U2 rU2 (bN2 )2 1+eU2 − μ − pN2 = r (1 − K ) − K − 1+ fU2 bU2 N2 V2 2 max{−pN2 , 1+ fU } = a, |Q21 | = (1+mN . fU2 )2 N2 2 V2 [2]

m1 (Q11 ) = j11 +

|Q12 | =

aV

− (1+eU2 )2 + 2

lU2 1+eU2

− μ − pN2 ,

To calculate m1 (Q22 ), the absolute value of the off-diagonal elements is added to the diagonal one in each column of Q22 and the maximum of the two sums is taken. Thus,



m1 (Q22 ) = max r 1 − + =

U2 K



−

N2 V2 rU2 aV2 (bN2 )2 mU2 − − + − ν + − K 1 + fU2 V2 N2 (1 + eU2 )2 1 + fU2

N V lU2 lV2 mU2 aU2 ; − μ − pN2 + −ν + 2 − 2 + 2 1 + fU2 V2 N2 1 + eU2 (1 + eU2 ) 1 + eU2



N V2 U2 mU2 − 2 + − ν + max r 1 − V2 N2 1 + fU2 K



−

rU2 (bN2 )2 − K 1 + fU2

aV2 lU2 lV2 aU2 − + ; − μ − pN2 + 1 + eU2 (1 + eU2 )2 (1 + eU2 )2 1 + eU2 =



N V2 mU2 − 2 + − ν + b, V2 N2 1 + fU2

where b = max{r (1 −



U2 K

 U2

)−

rU2 K

−

(bN2 )2 1+ fU2

aV

lV

− (1+eU2 )2 + (1+eU2 )2 ; 2 2

lU2 1+eU2

− μ − pN2 +

aU2 1+eU2 }.

rU2 aV2 (bN2 )2 lU2 − − − + − μ − pN2 + a, K K 1 + fU2 (1 + eU2 )2 1 + eU2 N V mU2 mN2 V2 g2 = 2 − 2 + −ν +b+ . V4 N2 1 + fU2 (1 + fU2 )2 N2 g1 = r 1 −



Therefore

(3.7)

Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

9

Rewriting (2.3) with I = 0 leads to

V2 lU2 = − pN2 − μ, V2 1 + eU2

(3.8)

N2 mU2 = − ν. N2 1 + fU2

(3.9)

Substituting (3.8) and (3.9) into (3.7) leads to





V2 V2 U2 rU2 aV2 (bN2 )2 +r 1− − − − + a ≤ − d, V2 K K 1 + fU2 V2 (1 + eU2 )2 V V mN2 V2 g2 = 2 + b + ≤ 2 − d, 2 V2 N (1 + fU2 ) 2 V2 g1 =

where d = min{−r (1 −

(bN2 )2

(3.10)

+ (1+eU2 )2 − a, − (1+ fU2 )2 N2 − b}. Therefore, m(Q ) ≤ 2 2 2 (3.10). Along each solution x(t, x0 ) to (2.3) such that x0 ∈ K and for t > T, we then have



1 t

t

0

U2 K

1 t

m(Q )ds ≤

)+



T

0

rU2 K

+

aV

mN

V

1+ fU2

m(Q )ds +

V2 V2

− d for t > T by (3.6) and

1 V (t ) t −T log −d t V (T ) t

which implies that m(Q ) ≤ −d/2 < 0. Hence, the following theorem is established . Theorem 1. The equilibrium E2 of the system (2.3) is globally asymptotically stable in  if m > fν , bN2 , 1+ fU2

m − f ν + eν > 0, β > δ +

(m − n )U2 1 + fU2

,

[2] d11

< 0,

[2] d22

< 0 and

[2] d33

ν

m− f ν

>

μ

l−eμ

, r (1 −

U2 K

)>

< 0.

3.3.3. Local and global stability of the steady state E3 For the equilibrium point E3 (U3 , 0, N3 , I3 ), the Jacobian matrix is given by

⎡

[3] j11

⎢ ⎢ 0 J3 = ⎢ ⎢ [m−(m−n)I3 ]N3 ⎣ (1+ fU3 )2

− 1+aUeU3 3

(c−1 )I3 ) − bU3 (1+ 1+ fU

[3] j22

0

0

[3] j33

0

0

m−n )U3 − ((1+ fU )2 3

[3]

where j11 = r − [3]

j33 =

2rU3 K

1+ fU3

⎥ ⎥ ⎥, (m−n )U3 N3 ⎥ −δ N3 − 1+ fU ⎦ 3 β − δ − (m1+−nfU)U3 3 0

bN3 (1+(c−1 )I3 ) , (1+ fU3 )2

[3]

j22 =

bian matrix J3 . From the second row this matrix has the eigenvalue are derived from the 3 × 3 Jacobian matrix given by

⎡

J3 =

⎢ [3] ⎣ f21 ⎡

=

[3] f11

[3] f12

[3] f13

[3] f22

[3] f23

0

[3] f33

[3] f31

3

⎢ [m−(m−n)I3 ]N3 ⎣ (1+ fU3 )2

[3] j33

[3]

[3]

[3]

⎤

bU3 N3 1+ fU3

−n )U3 N3 ⎥ −δ N3 − (m1+ ⎦, fU 3

β −δ−

0

3

[3]

− μ − pN3 (1 − I3 ), and the remaining eigenvalues

⎥ ⎦

(c−1 )I3 ) − bU3 (1+ 1+ fU

[3] j11

lU3 1+eU3

⎤

m−n )U3 − ((1+ fU )2

[3]

⎤

bU3 N3 1+ fU3

lU3 − μ − pN3 (1 − I3 ) and 1 + eU3 − δ I3 − ν . The local stability of the equilibrium at E3 is established from the eigenvalues of the Jaco-

−

(m−(m−n )I3 )U3

3

[3]

[3]

(m−n )U3 1+ fU3

[3]

[3]

[3]

[3]

[3]

[3]

[3]

[3]

[3]

where f11 = j11 , f12 = j13 < 0, f13 = j14 > 0, f21 = j31 > 0, f22 = j33 , f23 = j34 < 0, f31 = j41 < 0, f33 = j44 . The characteristic polynomial of the Jacobian of the system evaluated at this point is given by

λ3 + h1 λ2 + h2 λ + h3 = 0,

(3.11)

where [3] [3] [3] h1 = −( f11 + f22 + f33 ) = −Tr(J3 ), [3] [3] [3] [3] [3] [3] [3] [3] [3] [3] h2 = f11 f22 + f11 f33 + f22 f33 − f12 f21 − f13 f31 , [3] [3] [3] [3] [3] [3] [3] [3] [3] [3] [3] [3] h3 = − f12 f23 f31 + f13 f22 f31 + f12 f21 f33 − f11 f22 f33 = −Det J3 .

[3]

[3]

[3]

Now, if f11 < 0, f22 < 0, f33 < 0, then h1 > 0, h2 > 0 and h3 > 0. Furthermore, [3] [3] [3] 2[3] [3] [3] [3] [3] [3] 2[3] [3] [3] [3] [3] [3] [3] h1 h2 − h3 = f11 f12 f21 − f11 f22 + f12 f21 f22 − f11 f22 + f11 f13 f31 + f12 f23 f31

Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC 10

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23 2[3] [3] [3] [3] [3] 2[3] [3] [3] [3] [3] [3] 2[3] [3] 2[3] − f11 f33 − 2 f11 f22 f33 − f22 f33 + f13 f31 f33 − f11 f33 − f22 f33

> 0, if

[3] [3] [3]

2[3] [3]

[3] [3] [3]

[3] 2[3]

[3] [3] [3]

2[3] [3]

[3] [3] [3]

2[3] [3]

[3] [3] [3]

[3] 2[3]

f11 f12 f21 − f11 f22 + f12 f21 f22 − f11 f22 + f11 f13 f31 − f11 f33 − 2 f11 f22 f33 − f22 f33 + f13 f31 f33 − f11 f33 −

[3] 2[3]

[3] [3] [3]

f22 f33 > − f12 f23 f31 . From the Routh–Hurwitz criterion, the characteristic equation (3.11) has negative real parts. Thus, the following result can be stated. [3]

[3]

[3]

Lemma 3. Then the equilibrium point E3 is locally asymptotically stable if f11 < 0, f22 < 0, f33 < 0 and pN3 (1 − I3 ). The second compound matrix corresponding to J3 is given by

⎡ ⎢



−n )U3 N3 −δ N3 − (m1+ fU

0

[3] −n )U3 j11 + β − δ − (m 1+ fU

(c−1 )I3 ) − bU3 (1+ 1+ fU

(m−n )I3 (1+ fU3 )2

[m−(m−n )I3 ]N3 (1+ fU3 )2

[3] −n )U3 j33 + β − δ − (m 1+ fU

J3[2] = ⎣

P (U3 , N3 , I3 ) = diag 1,

⎡

0

⎢ ⎣

PF = ⎢0

N3 I3



0 and

0 N3 N3

⎡

0

−

⎢

N3 N3

0

I3 I3

I3 I3

⎢



N3 N3

0

⎥ ⎦.

−

I3 I3

⎤

1

0

0

P −1 = ⎣0

I3 N3

0⎦

0

0

I3 N3

⎢

I3 I3

−

⎤

⎡

⎥ ⎥, ⎦ 

0

0 N3 N3

0

⎡

3

⎤

0

  N3 I3

−

⎥ ⎦.

3



0

PF P −1 = ⎣0



3

N3 N3 , . I3 I3

0

Now

bU3 N3 − 1+ fU 3

The matrix function P is set as

Then

⎤

[3] [3] j11 + j33

3

⎥

⎤

[3] [3] j11 + j33

−n )U3 I3 −δ I3 − (m1+ fU

0

[3] −n )U3 j11 + β − δ − (m 1+ fU

(c−1 )I3 ) − bU3 (1+ 1+ fU

(m−n )N3 (1+ fU3 )2

[m−(m−n )I3 ]N3 (1+ fU3 )2

[3] −n )U3 j33 + β − δ − (m 1+ fU

PF J3[2] P −1 = ⎣

lU3 <μ+ 1 + eU3

bU3 I3 − 1+ fU

3

3

3

⎥ ⎦.

3

3

 [2]

The matrix Q = PF P −1 + PF J3 P −1 can be written in block form:



Q=

Q11 Q21

[3]



Q12 , Q22 [3]

with Q11 = j11 + j33 = r (1 −

U3 K

 −n )U3 I3 Q12 = −δ I3 − (m1+ fU 3

 Q22 =

rU3 K

)−

−

bN3 (1+(c−1 )I3 ) (1+ fU3 )2





bU3 I3 − 1+ , fU

Q21 =

3

[3] −n )U3 j11 + β − δ − (m + 1+ fU 3

[m−(m−n )I3 ]N3 (1+ fU3 )2

N3 N3

+

−

I3 I3

(m−(m−n )I3 )U3 1+ fU3

0 (m−n )N3 (1+ fU3 )2



− δ I3 − ν,

,



(c−1 )I3 ) − bU3 (1+ 1+ fU 3

[3] −n )U3 j33 + β − δ − (m + 1+ fU 3

N3 N3

−

I3 I3

.

A norm in R3 is defined as |u, n, i| = max{|u|, |n| + |i|} for any vector (u, n, i ) ∈ R3 . Let m denote the Lozinski ˘i measure with respect to this norm. Then

m(Q ) ≤ sup{ p1 , p2 },

(3.12)

with p1 = m1 (Q11 ) + |Q12 | and p2 = |Q21 | + m1 (Q22 ). where |Q12 | and |Q21 | are matrix norms induced by the L1 vector norm, and m1 denotes the Lozinski ˘i measure with respect to the L1 norm. More specifically, Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23 bN3 (1+(c−1 )I3 ) (m−(m−n )I3 )U3 + − δ I3 − ν, 1+ fU3 (1+ fU3 )2 (m−n )N3 = a, |Q21 | = (1+ fU )2 . Thus, 3

Now, as in Section 3.3.4, m1 (Q11 ) = r (1 − (m−n )U3 I3

|Q12 | = max{−δ I3 −

1+ fU3



m1 (Q22 ) = max r 1 − +

bU I , − 1+ 3fU3 } 3

U3 K



−

p2 =

−







−

rU3 bN3 (1 + (c − 1 )I3 ) − K (1 + fU3 )2

[m − (m − n )I3 ]N3 (m − (m − n )I3 )U3 bU3 (1 + (c − 1 )I3 ) ; − δ I3 − ν + 1 + fU3 1 + fU3 (1 + fU3 )2



I N3 (m − n )U3 − 3 +β −δ− + b, N3 I3 1 + fU3 U3 K

)−

rU3 K

−

Therefore



rU3 K

[m − (m − n )I3 ]N3 (m − (m − n )I3 )U3 (m − n )U3 N3 I3 ; − δ I3 − ν + β − δ − + − 2 1 + fU3 1 + fU3 N3 I3 (1 + fU3 )

where b = max{r (1 −

p1 = r 1 −

)−

I N3 U3 (m − n )U3 − 3 +β −δ− + max r 1 − N3 I3 1 + fU3 K +

=

U3 K

rU3 bN3 (1 + (c − 1 )I3 ) (m − n )U3 N3 I3 − +β −δ− + − 2 K 1 + fU3 N3 I3 (1 + fU3 )

bU3 (1 + (c − 1 )I3 ) + 1 + fU3 =

11

U3 K



−

bN3 (1+(c−1 )I3 ) (1+ fU3 )2

+

[m−(m−n )I3 ]N3 (m−(m−n )I3 )U3 ; 1+ fU3 (1+ fU3 )2

− δ I3 − ν +

bU3 (1+(c−1 )I3 ) }. 1+ fU3

(m − (m − n )I3 )U3 rU3 bN3 (1 + (c − 1 )I3 ) − + − δ I3 − ν + a, K 1 + fU3 (1 + fU3 )2

I N3 (m − n )U3 (m − n )N3 − 3 +β −δ− +b+ . N3 I3 1 + fU3 (1 + fU3 )2

(3.13)

Rewriting (2.3) with V = 0 leads to

N3 [m(1 − I3 ) + nI3 ]U3 = − δ I3 − ν, N3 1 + fU3 I3 = (1 − I3 ) I3

(3.14)

  (n − m )U3 β+ −δ .

(3.15)

1 + fU3

Substituting (3.14) and (3.15) into (3.13) gives



N3 U3 +r 1− N3 K



N3 rU3 bN3 (1 + (c − 1 )I3 ) − + a ≤ − d, K N3 (1 + fU3 )2 N I3 (m − n )N3 N3 p2 = 3 + +b+ ≤ − d, N3 1 − I3 (1 + fU3 )2 N3 p1 =

where d = min{−r (1 −

−

I bN3 (1+(c−1 )I3 ) ) − a, − 1−I3 (1+ fU3 )2 3

(m−n )N

− (1+ fU )32 − b}. Therefore, m(Q ) ≤ 3 (3.16). Along each solution x(t, x0 ) to (2.3) such that x0 ∈ K and for t > T, the result is

1 t

 0

t

m(Q )ds ≤

U3 K

1 t

)+

 0

T

rU3 K

+

m(Q )ds +

(3.16) N3 N3

− d for t > T by (3.12) and

1 N (t ) t −T log −d , t N (T ) t

which implies that m(Q ) ≤ −d/2 < 0. Hence, the following theorem is established. Theorem 2. The equilibrium point E3 of the system (2.3) is globally asymptotically stable in  if m > n, β > δ , (m − n ) > f (β − δ ), m(β − δ ) > ν (m − n ), C < 0 and G < 0.

3.3.4. Local and global stability of the interior steady state E4 The stability properties of the interior steady state in the presence of the disease in predator are studied. The local stability of E4 (U4 , V4 , N4 , I4 ) is established using the Routh–Hurwitz criterion. The Jacobian matrix of the system (2.3) at E4 is Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC 12

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

⎡

[4] j11

⎢ ⎢ (1+lVeU4 4 )2 J4 = ⎢ ⎢ (m−(m−n)I4 )N4 ⎣ (1+ fU4 )2

− 1+aUeU4 4

(c−1 )I4 ) − bU4 (1+ 1+ fU

[4] j22

−pV4 (1 − I4 )

m−n )U4 − ((1+ fU4 )2

where [4] j11 = r −

2rU4 K [4] j12

−

4

0

(m−(m−n )I4 )U4

0

0

1+ fU4

bN4 (1+(c−1 )I4 ) aV − (1+eU4 )2 (1+ fU4 )2 4 [4] [4] [4] 0, j13 < 0, j14 > 0, j21

Note that < of the linearized system is given by

− δ I4 − ν

[4]

and j22 = [4]

> 0, j23 <

⎤

bU4 N4 1+ fU4

⎥ ⎥ ⎥, (m−(m−n )U4 )N4 ⎥ −δ N4 − ⎦ 1+ fU4 β − δ − (m1+−nfU)U4 4 pN4V4

lU4 1+eU4 − μ − pN4 (1 − I4 ). [4] [4] [4] 0, j24 > 0, j31 > 0, j34 <

[4]

0 and j41 < 0. The characteristic polynomial

λ4 + m1 λ3 + m2 λ2 + m3 λ + m4 = 0,

(3.17)

where [4] [4] [4] [4] m1 = −( j11 + j22 + j33 + j44 ) = −Tr(J4 ), [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] m2 = j11 j22 + j11 j33 + j11 j44 + j22 j33 + j22 j44 + j33 j44 − j12 j21 − j13 j31 − j14 j41 , [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] m3 = j13 j22 j31 + j12 j21 j33 − j11 j22 j33 − j22 j33 j44 + j12 j21 j44 − j12 j24 j41 − j12 j23 j31 [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] − j11 j22 j44 − j13 j34 j41 + j13 j31 j44 + j14 j22 j41 + j14 j33 j41 − j14 j21 j42 − j11 j33 j44 , [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] m4 = − j12 j23 j34 j41 + j12 j24 j33 j41 + j12 j23 j31 j44 + j13 j34 j41 j22 − j13 j31 j22 j44 [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] − j14 j22 j33 j41 − j12 j21 j33 j44 + j11 j22 j33 j44 .

[4]

[4]

[4]

[4]

When j11 < 0, j22 < 0, j33 < 0 and j44 < 0, then m1 > 0 and m2 > 0. Furthermore, m3 > 0 if [4] [4] [4]

[4] [4] [4]

[4] [4] [4]

[4] [4] [4]

[4] [4] [4]

[4] [4] [4]

[4] [4] [4]

j13 j22 j31 + j12 j21 j33 − j11 j22 j33 − j22 j33 j44 + j12 j21 j44 − j11 j22 j44 − j13 j34 j41

[4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] + j13 j31 j44 + j14 j22 j41 + j14 j33 j41 − j11 j33 j44 > j12 j24 j41 + j12 j23 j31 , and m4 > 0 if [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] − j14 j22 j33 j41 − j12 j21 j33 j44 + j11 j22 j33 j44 − j13 j31 j22 j44 + j13 j34 j41 j22 > j12 j23 j34 j41 [4] [4] [4] [4] j12 j23 j31 j44 .

Let

[4] [4] [4] [4]

− j12 j24 j33 j41 −



[4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] q = ( j11 + j22 + j33 + j44 ) − j12 j21 − j13 j31 + j22 j33 − j14 j41 + j22 j44 + j33 j44



[4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] + j11 ( j22 + j33 + j44 )) j11 j22 j33 − j14 j22 j41 − j14 j33 j41 + j11 j22 j44 [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] + j11 j33 j44 + j22 j33 j44 − j13 ( j22 j31 − j34 j41 + j31 j44 )

 

[4] [4] [4] [4] [4] [4] [4] [4] + j12 ( j23 j31 + j24 j41 − j21 ( j33 j44 )) −

[4] [4] [4] [4] [4] [4] [4] [4] [4] j11 j22 j33 − j14 j22 j41 − j14 j33 j41

[4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] + j11 j22 j44 + j11 j33 j44 + j22 j33 j44 − j13 ( j22 j31 − j34 j41 + j31 j44 ) [4] [4] [4] [4] [4] [4] [4] [4] − j12 (− j23 j31 − j24 j41 + j21 ( j33 + j44 ))

2



[4] [4] [4] [4] 2 [4] [4] [4] [4] − ( j11 + j22 + j33 + j44 ) − j14 j22 + j33 j41



[4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] [4] + j22 ( j11 j33 j44 + j13 ( j34 j41 − j31 j44 )) + j12 ( j24 j33 j41 − j21 j33 j44 + j23 (− j34 j41 + j31 j44 )) .

From the Routh–Hurwitz criteria for a 4 × 4 Jacobian matrix, all the real parts of roots of (3.17) are negative if and only if the following conditions are satisfied,

m1 > 0, m3 > 0, m4 > 0 Now this is true if

[4] j11

< 0,

[4] j22

q = m1 m2 m3 − m23 − m21 m4 > 0.

and < 0,

[4] j33

< 0 and

[4] j44

(3.18)

< 0. Thus, the following result can be stated.

Lemma 4. The equilibrium point E4 is locally asymptotically stable if (a) r < (b) (c) (d)

2rU4 K

+

bN4 (1+(c−1 )I4 ) (1+ fU4 )2

aV

+ (1+eU4 )2 , 4 < μ + pN4 (1 − I4 ),

lU4 1+eU4 (m−(m−n )I4 )U4 < δ I4 1+ fU4 (m−n )U4 β < δ + 1+ fU . 4

+ ν,

The interior equilibrium is now shown to be globally asymptotically stable. ∗

∗

∗

bN [(1−I )+cI ] r Theorem 3. The interior equilibrium is globally asymptotically stable if M = (1+eUaV )(1+eU ∗ ) + (1+ fU )(1+ fU ∗ ) < K , j11 < 0, j22 < 0, lU4 [4] [4] j33 < 0, j44 < 0, μ < 1+eU , (nβ − mδ ) + (n − m )ν < 0, U4 < K, I4 > β (1 + cI4 )(1 − eU4 ), in which case V4 > 0. ∗

[4]

[4]

4

Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

13

Proof. A candidate Lyapunov function that has been used by Hsu [33], Dubey and Upadhyay [34], and others is chosen. The following function is considered,

  U    V  α1 U − U ∗ − U ∗ ln ∗ + α2 V − V ∗ − V ∗ ln ∗ U V   N    I  + α3 N − N ∗ − N ∗ ln + α4 I − I∗ − I∗ ln ∗ . ∗

Z (U, V, N, I ) =

N

(3.19)

I

Define Z1 (U ) = U − U ∗ − U ∗ ln( UU∗ ), Z2 (V ) = V − V ∗ − V ∗ ln( VV∗ ), Z3 (N ) = N − N ∗ − N ∗ ln( NN∗ ) and Z4 (I ) = I − I∗ − I∗ ln( II∗ ). Z can be rewritten as Z (U, V, N, I ) = α1 Z1 (U ) + α2 Z2 (V ) + α3 Z3 (N ) + α4 Z4 (I ). Differentiating Z with respect to time t along the solutions of model (2.3) yields

dZ = dt

        U ∗ dU V ∗ dV N ∗ dN I∗ dI α1 1 − + α2 1 − + α3 1 − + α4 1 − . U

dt

V

dt

N

dt

I

(3.20)

dt

The linear approximations U − U ∗ ∼ = 1 + eU ∼ = 1 + fU, V − V ∗ ∼ = V, N − N ∗ ∼ = N and I − I∗ ∼ = I are used to compute dZ2 (V (t )) dZ3 (N (t )) , dt dt

dZ1 = dt



dZ4 (I (t )) dt

and

1−

  U∗

as follows:

r 1−

U

dZ1 (U (t )) , dt

U K



−



aV bN[(1 − I ) + cI] − U 1 + eU 1 + fU

r a(U − U ∗ )(V − V ∗ )(1 + eU ∗ ) − aV ∗ (U − U ∗ )2 = − (U − U ∗ )2 − K (1 + eU )(1 + eU ∗ ) − dZ2 = dt



b(N − N ∗ )[(1 − I ) + cI](1 + fU ∗ )(U − U ∗ ) − bN ∗ [(1 − I∗ ) + cI∗ ](U − U ∗ )2 , (1 + fU )(1 + fU ∗ )

1−

V∗ V

 lU



1 + eU

= (V − V ∗ )

 lU

− p( 1 − I ) N − μ V

−

lU ∗ lU ∗ + − p( 1 − I ) N − p ( 1 − I ∗ ) N ∗ + p ( 1 − I ∗ ) N ∗ − μ 1 + eU ∗ 1 + eU ∗



1 + eU l[U ∗ − (1 + eU ∗ ) =− (U − U ∗ )(V − V ∗ ) − p(V − V ∗ )[(1 − I )N − (1 − I∗ )N∗ , (1 + eU )(1 + eU ∗ ) dZ3 = dt



1−

N∗ N

 [m(1 − I ) + nI]U 1 + fU



− δI − ν N

 [m(1 − I ) + nI](U − U ∗ )(1 + fU ∗ ) − [m(1 − I∗ ) + nI∗ ]U ∗ (U − U ∗ )  = (N − N ) − δ (I − I ∗ ) ∗ (1 + fU )(1 + fU ) [m(1 − I ) + nI](1 + fU ∗ ) − [m(1 − I∗ ) + nI∗ ]U ∗ = (U − U ∗ )(N − N∗ ) − δ (I − I∗ )(N − N∗ ) (1 + fU )(1 + fU ∗ ) [m(1 − I∗ ) + nI∗ ]U ∗ − [m(1 − I ) + nI](1 + fU ∗ ) =− (U − U ∗ )(N − N∗ ) − δ (I − I∗ )(N − N∗ ) (1 + fU )(1 + fU ∗ ) ∗

and

dZ4 = dt



1−



I∗ (1 − I ) I

  (m − n )U β− −δ I 1 + fU

 (m − n )U (m − n )U ∗ (m − n )U ∗  = (I − I∗ )(1 − I ) β − δ − + − ∗ ∗ 1 + fU

1 + fU

1 + fU

(m − n )(U − U ∗ )(1 + fU ∗ )(I − I∗ )(1 − I ) + (m − n )U ∗ (U − U ∗ )(I − I∗ )(1 − I ) =− (1 + fU )(1 + fU ∗ ) ∗ ∗ (1 + fU ) − U = − (m − n ) (U − U ∗ )(I − I∗ )(1 − I ). (1 + fU )(1 + fU ∗ ) Now,

dZ = dt

 r a(U − U ∗ )(V − V ∗ )(1 + eU ∗ ) − aV ∗ (U − U ∗ )2 α1 − (U − U ∗ )2 − K (1 + eU )(1 + eU ∗ ) b(N − N ∗ )[(1 − I ) + cI](1 + fU ∗ )(U − U ∗ ) − bN ∗ [(1 − I∗ ) + cI∗ ](U − U ∗ )2 − (1 + fU )(1 + fU ∗ )



+ α2 −



l[U ∗ − (1 + eU ∗ ) (U − U ∗ )(V − V ∗ ) − p(V − V ∗ )[(1 − I )N − (1 − I∗ )N∗ (1 + eU )(1 + eU ∗ )



Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC 14

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

 + α3 −

[m(1 − I∗ ) + nI∗ ]U ∗ − [m(1 − I ) + nI](1 + fU ∗ ) (U − U ∗ )(N − N∗ ) (1 + fU )(1 + fU ∗ )





 (1 + fU ∗ ) − U ∗ (U − U ∗ )(I − I∗ )(1 − I ) ∗ (1 + fU )(1 + fU )  r  ∗ ∗ aV bN [(1 − I∗ ) + cI∗ ] = α1 − + + (U − U ∗ )2 K (1 + eU )(1 + eU ∗ ) (1 + fU )(1 + fU ∗ )   a(1 + eU ∗ ) l[U ∗ − (1 + eU ∗ ) + −α1 − α (U − U ∗ )(V − V ∗ ) 2 (1 + eU )(1 + eU ∗ ) (1 + eU )(1 + eU ∗ )   b[(1 − I ) + cI](1 + fU ∗ ) [m(1 − I∗ ) + nI∗ ]U ∗ − [m(1 − I ) + nI](1 + fU ∗ ) + −α1 − α3 (U − U ∗ )(N − N∗ ) ∗ ∗ (1 + fU )(1 + fU ) (1 + fU )(1 + fU ) − α2 p(V − V ∗ )[(1 − I )N − (1 − I∗ )N ∗ ] − α3 δ (I − I∗ )(N − N ∗ )   (1 + fU ∗ ) − U ∗ − α4 (m − n ) (1 − I ) (U − U ∗ )(I − I∗ ). ∗ (1 + fU )(1 + fU ) − δ (I − I∗ )(N − N ∗ ) + α4 −(m − n )

∗

∗

∗

∗

(1+eU )−U (1+eU ) a (1+eU ) Let α1 = 1 and α4 = 1, then α2 (l[1+ = α1 (1+aeU eU )(1+eU ∗ ) )(1+eU ∗ ) , that is, α2 = l[(1+eU ∗ )−U ∗ . Also

α3

b[(1 − I ) + cI](1 + fU ∗ ) [m(1 − I∗ ) + nI∗ ]U ∗ − [m(1 − I ) + nI](1 + fU ∗ ) =− , ∗ (1 + fU )(1 + fU ) (1 + fU )(1 + fU ∗ )

that is, α3 =

b[(1−I )+cI](1+ fU ∗ ) . [m (1−I )+nI](1+ fU ∗ )−[m (1−I∗ )+nI∗ ]U ∗

Hence,

 r



dZ aV ∗ bN ∗ [(1 − I∗ ) + cI∗ ] = − + + (U − U ∗ )2 dt K (1 + eU )(1 + eU ∗ ) (1 + fU )(1 + fU ∗ ) a(1 + eU ∗ ) − p(V − V ∗ )[(1 − I )N − (1 − I∗ )N ∗ ] l[(1 + eU ∗ ) − U ∗ b[(1 − I ) + cI](1 + fU ∗ ) − δ (I − I∗ )(N − N∗ ) [m(1 − I ) + nI](1 + fU ∗ ) − [m(1 − I∗ ) + nI∗ ]U ∗   (1 + fU ∗ ) − U ∗ − (m − n ) (1 − I ) (U − U ∗ )(I − I∗ ). ∗ (1 + fU )(1 + fU ) The coefficient of (U − U ∗ )2 is strictly negative if

−

r aV ∗ bN ∗ [(1 − I∗ ) + cI∗ ] + + < 0, K (1 + eU )(1 + eU ∗ ) (1 + fU )(1 + fU ∗ )

that is,

M=

aV ∗ bN ∗ [(1 − I∗ ) + cI∗ ] r + < . ∗ (1 + eU )(1 + eU ) (1 + fU )(1 + fU ∗ ) K

If α 2 > 0 and α 3 > 0, then the function Z is negative definite.



3.4. Persistence of the system A system is defined to be uniformly persistence if the minimum of each component u(t) of a positive solution is always greater than some positive constant, i.e. lim inf u(t ) >  > 0 [35]. The persistence of a system implies that all the species t→∞

continue to exist and none of them will become extinct [34]. From Proposition 3 the system (2.3) has been proved to be bounded. Now, the system is shown to be persistent. To do so, all the boundary equilibria are shown to be repellers. Theorem 4. If the following conditions hold (i) β > δ , lU3 (ii) 1+eU > μ + pN3 (1 − I3 ), 3

then the system is persistent. [1] Proof. From the Jacobian matrix associated with E1 (U1 , 0, N1 , 0) the following eigenvalue λ1 = β − δ − ν + nmν = β − δ +

ν

was found. But equilibrium point E3 exists if m > n and β > δ , making λ1 > 0. Hence, existence of E3 implies that E1 is unstable. (m−n )U For the equilibrium point E2 (U2 , V2 , N2 , 0), from the Jacobian matrix J2 , there is an eigenvalue β − δ − 1+ fU 2 = m (n − m )

β − δ − (m − n ) mν . Using the same argument as above, E2 is unstable if E3 exists.

[1]

2

Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

15 [3]

For the equilibrium point E3 (U3 , 0, N3 , I3 ), from the Jacobian matrix J3 , an eigenvalue j22 =

lU3 1+eU3

− μ − pN3 (1 − I3 )

exists. If > μ + pN3 (1 − I3 ) then the eigenvalue is positive, and E3 is unstable. Thus, all the boundary equilibria of system (2.3) are repellers if the conditions stated in the theorem hold. lU3 1+eU3



3.5. The effect of super-predator removal on the eco-epidemiological system In the absence of the super-predator, the system (2.2) reduces to a prey–predator sub-system given by







dU U aV = U r 1− − , dt K 1 + eU   dV lU =V −μ . dt 1 + eU

(3.21)

This subsystem has three biologically feasible equilibria namely: (i) E¯0 (0, 0 ), (ii) E¯1 (K, 0 ) and (iii) E¯2 (U2 , V2 ) where U2 = and V2 =

rl (Kl −μ−eK μ ) ). aK (eμ−l )2

The first two equilibria always exist, and E¯2 exists if l > eμ and l − eμ >

of the system (3.21) is given by



J (U, V ) =

r (1 −

2U K

)− lV

μ K

μ

l−eμ

or U2 < K. The Jacobian



aV

− 1+aUeU . lU −μ 1+eU

(1+eU )2

(1+eU )2

(3.22)

(a) E¯0 is always unstable with eigenvalues r and −μ. It is a saddle point whose stable manifold is the V-axis. ¯ (b) E¯1 has eigenvalues −r and 1+KleK − μ. Thus, E¯1 is unstable in the V-axis if 1+KleK > μ, that is, l − eμ > μ K and E2 exists. μ ¯ It is locally asymptotically stable if l − eμ < K and E2 does not exist. Next, the global stability of E¯1 is proved. In order to do so, Theorem 7 and Lemma 6 that are essential in finding a candidate Lyapunov function and proving the global stability are stated. Define L : {(U, V ) ∈  : U > 0} → R by L(U, V ) = 12 V 2 . The time derivative of L computed along solutions of (3.21) is

L (U, V ) = V · V  = V 2

 ( l − eμ ) −

μ U

U (1 + eU )



.

It follows that L (U, V) ≤ 0 if U = K, that is, l − eμ − μ K < 0. Hence, L is a Lyapunov function on . Substituting V = 0 in the first equation of (3.21) leads to dU = Ur (1 − UK ), which shows that U → K as t → ∞. Therefore, it follows from dt the LaSalle’s Invariance Principle, that every solution of the equations in the model (3.21), with initial conditions in , approaches E¯1 as t → ∞. Hence, it follows that the point E¯1 is globally asymptotically stable in  if l − eμ − μ K < 0.

(c) E¯2 has eigenvalues λ1 =



−a1 + a21 −4a0 a2 2a0

−a1 − a21 −4a0 a2 and λ2 = , the 2a0 lrμ − eKlrμ + erμ2 + e2 Krμ2 and

roots of the quadratic equation a0 λ2 + a1 λ + a2 =

0, where a0 = Kl (l − eμ ), a1 = a2 = Kl 2 rμ − lrμ2 − 2eKlrμ2 + erμ3 + e2 Krμ3 . Now, a0 > 0, since it is derived from a condition for existence of E¯2 . a1 = rμ[(l + eμ ) − eK (l − eμ )]. a1 > 0 if μ

eK. a2 = (l − eμ )r K [(l − eμ ) −

μ K

]. Thus, a2 > 0 if l − eμ >

μ K

l+eμ l−eμ

>

. The following result can be stated.

Lemma 5. E¯2 is locally asymptotically stable when it exists if a1 > 0, a2 > 0, that is,

l+eμ l−eμ

> eK and l − eμ > μ K.

Furthermore,

Tr(J (U2 , V2 )) = −μ − −

rμ elrμ(Kl − μ − eK μ ) lμ + eμ 2 + K ( l − eμ ) K (l − eμ )3 (1 + l−e ) ( l − e μ )( 1+ μ

lr (Kl − μ − eK μ ) K ( l − eμ )

Det(J (U2 , V2 )) = −r μ −



2

1+

eμ l−eμ



 +r 1+

μ

K ( l − eμ )



eμ l−eμ

)

.

el 2 r μ2 (l − eμ − μ ) lrμ(l − eμ − 3Kμ ) 2 r μ2 lrμ K − + + e μ eμ K ( l − eμ ) (l − eμ )3 (1 + l−eμ )2 (l − eμ )2 (1 + l−e ) ( l − e μ )( 1+ μ

eμ l−eμ

)

.

The type of interior equilibrium point depends on the sign of

(U¯ 2 , V¯2 ) = Tr2 (J (U2 , V2 )) − 4Det(J (U2 , V2 )) =

1 (rμ(rμ(l + eμ )2 − K 2 (l − eμ )2 (4l 2 − 4el μ − e2 rμ ) K 2 l 2 ( l − eμ ) 2 + 2K μ(l − eμ )(2l 2 − e2 rμ − el (r + 2μ )))).

If (U¯ 2 , V¯2 ) ≥ 0 then the interior equilibrium is a stable node, otherwise it is a spiral point. The existence condition of E¯2 renders l > eμ. Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC 16

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

Let δ (x ) = d2 x2 + d1 x + d0 where d2 = −K 2 (4l 2 − 4el μ − e2 rμ ), d1 = 2K μ(2l 2 − e2 rμ − el (r + 2μ )) and d0 = rμ(l + eμ )2 . d0 > 0. The discriminant

d12 − 4d0 d2 = 4K 2 rμ(l + eμ )2 (4l 2 − 4el μ − e2 rμ ) + 4K 2 μ2 (2l 2 − e2 rμ − el (r + 2μ ))2 > 0 if 4l (l − eμ ) > e2 rμ. Now, the global stability of E¯2 is discussed whenever it exists, using the Bendixson–Dulac Criterion [25] and the Poincar e´ –Bendixson Theorem [26]. Proposition 5. E¯2 is globally asymptotically stable if 2μ > lK. Proof. Consider the system (3.21). Let (U, V ) = in the interior of the U − V plane. Then

1 UV

, h1 (U, V ) = Ur (1 −

U K

)−

aUV 1+eU ,

h2 (U, V ) =

lUV 1+eU

− μV . Now (U, V) > 0

∂ ( h1 ) ∂ ( h2 ) + ∂U ∂V     ∂  l U a ∂ r μ = 1− − + − ∂U V K 1 + eU ∂ V 1 + eU U aeKV − r (1 + eU )2 = . KV (1 + eU )2

(U, V ) =

Now  is not identically zero in the positive quadrant of the U − V plane. (U, V) < 0 and does not change sign if

aeKV − r (1 + eU )2 = aeKV − r − 2reU − e2U 2 = −r − (2rU − aKV )e − e2U 2 < 0, that is, if 2rU − aKV > 0. This means that at E2 , we have

2rU2 − aKV2 =

rl (Kl − μ − eK μ ) 2r μ − l − eμ ( l − eμ ) 2

r (l − eμ )(2μ − lK ) + μrl ( l − eμ ) 2 > 0, =

¯ . Since E¯2 is locally asymptotically stable, the if and only if 2μ > lK. Thus, subsystem (3.21) does not have a limit cycle in  Poincar e´ –Bendixson Theorem and Proposition 3 imply that the equilibrium point E¯2 is globally asymptotically stable in the ¯.  region  4. Numerical analysis It was estimated that 1684 lions lived in the KNP in 2005 and 2006 [36]. The lion population in KNP appears to be stable at around 1700 since 2005. The KNP has a current lion population of approximately 1700. Keet et al. [37] estimated that about 500 of the 1700 lions reside in areas with buffalo whose BTB prevalence is high [20]. The initial population densities for each species are chosen to be positive at U (0 ) = 133, 0 0 0 [36], V (0 ) = 412 [38], N (0 ) = 1700 [36], W (0 ) = 1200 and B(0 ) = 500. Some parameter values which were not readily available were computed on the basis of the steady states and stability conditions, and these include prey conversion rates into predator biomass. 4.1. Parameter estimates The literature on the dynamics of the three species in KNP as well as studies from other parks with similar environmental conditions (e.g. the Serengeti National Park SNP) provide the following information: (i) Carrying capacity for impala The impala population was estimated to lie between 132,300 and 176,400 between the years 2010 and 2011, [36]. In this study, the carrying capacity was estimated to be 20 0,0 0 0. (ii) In a study on impala in KNP, Fairall [39] deduced that fecundity lies around 95% in mature females and is drastically low in two-year old females. A sample of 100 females studied produced 244 lambs during their lifetime. A female was observed to produce offspring for a period of 10 years up to the age of 12 years. At an early reproducing age of two years, approximately 45% of the mortality had taken place. In this study, the growth rate of the impala was estimated to be 0.01342. (iii) Mortality rates The life-expectancy of cheetah in protected areas is estimated at 18 years [41]. In this study the mortality rate of 1 cheetah was estimated as the reciprocal of the life-expectancy. Hence, μ = 18 = 0.0556. Once the juvenile cheetah emerge from the lair, mortality rate due to predation by lions is estimated at 50% [40]. In this study the cheetah mortality rate by the lions was estimated as p = 0.5. The lifespan for lions in KNP is about 15 years [23]. The mortality 1 rate of lions was estimated as the reciprocal of the life-expectancy. Hence, ν = 15 = 0.0667. Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

ARTICLE IN PRESS

JID: AMC

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

17

Table 2 Parameter description and values. Symbol

Description

Value

Reference

r K a b c e f l m

Intrinsic growth rate of the impala Environmental carrying capacity of the impala Capture rate of the impala by the cheetah Capture rate of the impala by the lion Efficiency of infected lion to capture impala Handling time of impala by cheetah Handling time of impala by lion Impala biomass conversion rate into new cheetah Impala biomass conversion rate into new susceptible lion Impala biomass conversion rate into new infected lion Mortality rate of cheetah by lion Disease standard incidence Natural mortality rate of lion Natural mortality rate of cheetah Disease-induced mortality rate of infected lion

0.01342 20 0,0 0 0 0.0 0 01787 0.0 0 0 09525 0.5 0.0 0 03 0.0 0 04 0.0 0 0 01705

[39] Estimate [15] [44] Estimate [15] Estimate Computed

0.0 0 0 027

Computed

0.0 0 0 0135 0.5 0.16 0.0667 0.0556 0.0767

Computed [40] [23] [23] [41] [23]

n p

β ν μ

δ

(iv) Predation rates In a study on feeding ecology of cheetah in south eastern KNP, Mills et al. [15] deduced that impala were the dominant prey killed out of nine prey species selected. Seven adult cheetahs were monitored in the same region between 1987 and 1990 [42]. Kill rates were observed for different coalitions. It was deduced that on average a cheetah consumes 1.4 kg of meat per day [42]. Out of the total weight of an adult impala, only 60% was estimated to be consumable [15]. The body masses of adult male and female impala were estimated as 54.4 kg and 40.9 kg [16,43], giving an average weight of 47.65 kg. Thus, an average cheetah eats 17.87 impala per year. The predation rate is considered as the number of prey killed/consumed per predator per year. In this study, the cheetah predation rate of the impala was considered as 0.0 0 01787 after re-scaling. The impala is the third dominant prey species consumed by the lion after zebra and wildebeest in the KNP [17]. Between 1986 and 1995, the following impala kill rates per lion per year were observed in the KNP: Pride females 13.0 (1986–1989), pride females 8.7 (1992–1995), territorial males 1.3 and non-territorial males 15.1 [44]. In this study it was estimated that the impala kill rate per lion per year to be the average 9.525. In this study, the lion predation rate of the impala was considered as 0.0 0 0 09525 after re-scaling. (v) Kill retention times Kill retention times are influenced by how large the group and prey are, densities of predators, and awareness of the presence of competing predators [15]. The cheetah kill retention time was 165 min [15]. In this study the 165 cheetah handling time considered as kill retention time per day was estimated to be e = 1440 = 0.1146 per day, or 165 e = 1440 = 0 . 0 0 03 per year. The lion is expected to have a longer retention time as it has no competitors. In this ×365 study the kill retention time was estimated to be 220 min, giving f = 0.1528 per day or f = 0.0 0 04 per year. (vi) Conversion rates The cheetah and lion biomass conversion rates l and m were estimated in relation with the bounds for existence of equilibrium points E2 and E3 . For E2 , l > eμ = 1.668 × 10−5 . For E3 , m > f ν = 2.668 × 10−5 . In this study, it was estimated that l = 1.705 × 10−5 and m = 2.7 × 10−5 . (vii) Disease standard incidence and disease-induced mortality rate “The initial frequency of diseased animals is defined as the proportion of animals in the population that contract the disease (with or without clinical symptoms) and are infected with the pathogen”, [23]. In a study by Keet (unpublished data 1999–2004), 16 animals tested positive and were observed for a period of five years. Of the 16, 7 out of 10 changed stage from infected to diseased, whilst one lioness succumbed to injuries though still infectious. Thus, 80% of infected lions became diseased within five years [22], leading to a deduction that 8% of lions became diseased within a period of six months. In this study, the initial frequency of diseased animals was estimated to be 0.16 per year. Furthermore, in the same study by Keet, the mortality from disease at the end of infectious period was 100%. 14 diseased lions were monitored to check the possibility of recovery. However, all succumbed to the disease in a period of six months [23]. In this study, the disease-induced mortality rate was estimated to be δ = 0.0767.

The parameters values used in the numerical analysis of model (2.2) and subsystem (3.21) are given in Table 2. Keeping the rest of the parameters at their fixed levels, some ecological and epidemiological parameters are varied, namely p, β and δ to try to understand their influence on the dynamics of system (2.2). The re-scaled species population in relation to its initial value, U(t)/U(0) for instance, was plotted as a function of time. Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

JID: AMC 18

ARTICLE IN PRESS

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

Fig. 1. (a) Scaled populations, (b) trajectory tending to E¯2 .

Fig. 2. (a) Scaled populations, (b) trajectory, when cheetah mortality by lion, p = 0 and all other parameter values remain unchanged.

4.2. The effect of super-predator removal on the system dynamics In the absence of the lion, impala and cheetah populations and phase diagrams are shown in Fig. 1. It can be observed that the positive equilibrium E¯2 is reached as time gets large. The impala population rises rapidly to reach a plateau below the carrying capacity, before it falls down and settle to equilibrium value above the initial level. The cheetah population rises to more than double the initial value. The trajectories indicate the global stability of the equilibrium point E¯2 .

4.3. The effect of cheetah mortality by lion on cheetah population density In this section the effect on the system of changes in the cheetah mortality by lion, p, is described. As the mortality rate increases in the wide range 0 ≤ p ≤ 0.5, the cheetah population becomes depressed and eventually goes extinct. For instance, in the absence of lion kills, p = 0, the cheetah population is stable. The population of the impala grows to below carrying capacity. However, the lion population becomes extinct in 50 years. Increasing p from 0 to 0.0 0 0 01, the cheetah population takes a slight setback and eventually stabilises at about 350. As p takes the value 0.001, the cheetah population drops to zero spontaneously. This shows the effect the cheetah mortality by lion has on the cheetah population. Thus, the added cheetah mortality by lions has the potential of eliminating the cheetah species from the system. There is only a quantitative change, qualitatively, the system’s behaviour remains the same. This behaviour is shown in Figs. 2–4. In all the three cases both lion populations decline and tend to zero. This occurs as a result of high disease incidence rate. Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

JID: AMC

ARTICLE IN PRESS

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

19

Fig. 3. (a) Scaled populations, (b) trajectory, when cheetah mortality by lion, p = 0.0 0 0 01 and all other parameter values remain unchanged.

Fig. 4. (a) Scaled populations, (b) trajectory, when cheetah mortality by lion, p = 0.001 and all other parameter values remain unchanged.

4.4. Effect of epidemiological parameters on the system dynamics In this section the effects of changes in the disease-related parameters, β and δ , on the system are described. The system changes are considered in the absence of cheetah kills by lions, p = 0, as the cheetah population stabilises and does not face extinction. All other parameter values remain constant. 4.4.1. Disease standard incidence in lions, β The disease standard incidence, β , is 0.16. At this value, the total lion population does not survive long enough, but declines continuously and becomes extinct in 50 years. When the magnitude of β is decreased from 0.16 to 0.08 the healthy lions survive but the infected lions becomes extinct in 150 years, twice the initial period. Thus, the disease is eliminated. The cheetah population survives. Fewer lions are becoming infected and as such the healthy lions survive and predation is enhanced. This behaviour is shown in Figs. 5 and 6. As such, a decrease in the disease incidence rate has a positive effect on the lion population, negative effect on the impala population, and has insignificant effect on the cheetah population. The effects remain the same when β is reduced further to β = 0.02. 4.4.2. Disease-induced death rate, δ The disease-induced mortality rate, δ , is 0.0767. When the value of δ is decreased from 0.0767 to 0.00767, the susceptible lion population are still eliminated from the system. The diseased lions survive longer. Since the rate of infection has not changed, more susceptible lions become diseased and are eliminated from the system. A further decrease in δ to 0.0 0 0767 is accompanied by only a marginal increase in the diseased lion population. However, the entire lion population still becomes extinct after 100 years. This behaviour is shown in Figs. 7 and 8. Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

JID: AMC 20

ARTICLE IN PRESS

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

Fig. 5. (a) Scaled populations, (b) trajectory, when disease incidence, β = 0.08 and p = 0 and all other parameter values remain unchanged.

Fig. 6. (a) Scaled populations, (b) trajectory, when disease incidence, β = 0.02 and p = 0 and all other parameter values remain unchanged.

Fig. 7. (a) Scaled populations, (b) trajectory, when disease induced mortality rate, δ = 0.00767 and p = 0 and all other parameter values remain unchanged.

Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

JID: AMC

ARTICLE IN PRESS

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

21

Fig. 8. (a) Scaled populations, (b) trajectory, when disease induced mortality rate, δ = 0.0 0 0767 and p = 0 and all other parameter values remain unchanged.

5. Discussion and conclusions In relation to cheetah mortality by the lions, p, the cheetah in KNP can only survive in the absence of lion killings, p = 0, or killings are kept at minimum levels p = 0.0 0 01. It is recommended that the two species be separated to ensure that the cheetah species does not become extinct. It appears that the current population level of the cheetah is too low to ensure growth beyond its initial value even with a situation in which there is no competition by lions. Thus, the absence of lion in the system allowed the cheetah population to rise significantly to more than double the initial value. The equilibrium point E¯2 is reached as shown in Fig. 1. The results indicate that cheetah thrives very well in the absence of large competition from the lion. These findings are in agreement with the studies by Linnell and Strand [45]. Removal of large predators from an ecosystem may allow small to medium predator density to improve. However, accompanied by this is the enhanced predation of prey populations [45–47]. The impala population rises to about 1.18 times the initial value, lower than the 1.4 times the initial value in the presence of the lion. This can be seen upon comparing Figs. 1–3. μ In the absence of disease, the three species coexist. The equilibrium point E2 exists if m > fν , m−ν f ν > l−e μ , that is, the growth ratio of super-predators exceed that of predators, r (1 −

U2 K

)>

bN2 , 1+ fU2

the growth function of prey exceeds their

predation rate by super-predators. Furthermore, it is locally asymptotically stable if β > δ +

(m − n )U2

[2]

[2]

, d11 < 0, d22 < 0 1 + fU2 [2] ˘ and d33 < 0. It is globally stable if the Lozinski i measure with respect to the L1 norm is such that m(Q ) ≤ −d/2 < 0. The value of the disease standard incidence, β , is high resulting in that the susceptible lions become extinct in about 25 years. At the same time the entire lion population becomes diseased. The prediction is that the entire lion population disappear in 50 years. When β = 0.08. the lion population persists but is not able to grow beyond its initial value. The value of the disease-induced mortality rate, δ , is high. All infected lions are dying early. At δ = 0.00767 the diseased lions survive longer but remain infectious, further spreading the disease to the healthy lions. The reduced mortality rate leads to an increase in diseased lions before they die. Cheetah and lion co-exist as long as disease infection rate assumes very low values. Cheetah population persists, but grows marginally with time. Furthermore, the presence of disease in lion reduces the pressure lion exerts on cheetah, resulting in cheetah able to survive longer, but only as long as added mortality is kept very low. This is shown in Figs. 2, 5 and 6. The equilibrium point E4 for the system (2.3) was found to be locally asymptotically stable, using Routh Hurwitz [4] [4] [4] [4] criterion, if the conditions j11 < 0, j22 < 0, j33 < 0 and j44 < 0 were satisfied. Furthermore, as in Theorem 3 global ∗

∗

∗

bN [(1−I )+cI ] r asymptotic stability was assured if M = (1+eUaV )(1+eU ∗ ) + (1+ fU )(1+ fU ∗ ) < K . A reduction in disease infection rate or disease induced death rate does not impart any increase in the cheetah population. As long as the disease is present and the disease incidence is high, the lion population remains suppressed and eventually becomes extinct. The conclusions from this study do not agree with Maas et al.’s [19] findings. They assessed the impact of feline immunodeficiency virus (FIV) and BTB co-infection in African lions in KNP. They used a multivariable logistic regression model for analysis. They found that BTB does not pose a serious conservation threat to KNP, and that a significant spatio-temporal increase of BTB may impact on lion health. The presence of disease in lions cannot be used as biological control on the cheetah population. However, the disease affects the lion species. At high levels of disease incidence, lion population becomes extinct rapidly. The continued survival ∗

Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

JID: AMC 22

ARTICLE IN PRESS

[m3Gsc;October 28, 2016;20:34]

W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

of the lions depends on the levels of the disease incidence and the disease-induced mortality rates. Without human intervention, the study showed the lions becoming extinct. Urgent treatment of disease is essential to reduce further infection and save the lion from possible extinction. The reduced infection rate is accompanied by a decline in prey population, as more healthy lions hunt. Acknowledgements The financial assistance of the National Research Foundation (NRF) (Grant no. 95380) towards this research is hereby acknowledged. Opinions expressed and conclusions arrived at, are those of the authors and are not necessarily to be attributed to the NRF. The authors appreciate the financial support provided by the Nelson Mandela Metropolitan University. Appendix A. Auxiliary and standard results of the theory of ordinary differential equations Theorem 5. (Bendixson–Dulac criterion) [25] Let h1 (u, v ), h2 (u, v ) and (u, v ) be C1 functions in a simply connected domain D ⊂ R2 such that

∂ ( h1 ) ∂ ( h2 ) + ∂u ∂v does not change sign in D and vanishes at most on set of measure zero. Then the system

du = h1 (u, v ) dt dv = h2 (u, v ) dt has no periodic orbits in D. Theorem 6. (Poincare´ –Bendixson) [26] Let the functions h1 and h2 have continuous first partial derivatives in a domain D of the uv-plane. Let D1 be a bounded sub-domain in D, and let R be the region that consists of D1 plus its boundary (all points of R are in D). Suppose that R contains no critical point of the system

du = h1 (u, v ) dt dv = h2 (u, v ) dt If there exists a constant t0 such that x = φ (t ), y = ψ (t ) is a solution of the system that exists and stays in R for all t ≥ t0 , then either x = φ (t ), y = ψ (t ) is a periodic solution (closed trajectory), or x = φ (t ), y = ψ (t ) spirals toward a closed trajectory as t → ∞. In either case, the system has a periodic solution in R. Theorem 7. [48] Let E be an open subset of Rn containing x0 . Suppose that f ∈ C1 (E) and that f(x0 ) = 0. Suppose further that there exists a real valued function V ∈ C1 (E) satisfying V (x0 ) = 0, and V(x) > 0 if x = x0 . Then (a) if V˙ (x ) < 0 for all x ∈ E, x0 is stable; (b) if V˙ (x ) < 0 for all x ∈ E ∼ {x0 }, x0 is asymptotically stable; (c) if V˙ (x ) > 0 for all x ∈ E ∼ {x0 }, x0 is unstable. Lemma 6. (Lasalle Invariance Principle) [33] Assume that V is a Lyapunov function of the dynamical system x˙ = f(x ) on E. Define S = {x ∈ G¯ ∩ E : V˙ (x ) = 0}. Let M be the largest invariant set in S. Then every bounded trajectory (for t ≥ 0) of x˙ = f(x ) that remains in G approaches the set M as t → ∞. References [1] R.A. Saenz, H.W. Hethcote, Competing species models with an infectious disease, J. Math. Biosci. Eng. 3 (1) (2006) 219–235. [2] M. Haque, A predator-prey model with disease in the predator species only, J. Nonlinear Anal. Real World Appl. 11 (2010) 2224–2236. [3] B. Mukhopadhyay, R. Bhattacharyya, Dynamics of a delay-diffusion prey-predator model with disease in the prey, J. Appl. Math. Comput. 17 (1–2) (2005) 361–377. [4] Y. Lu, D. Li, S. Liu, Modeling the hunting strategies of the predators in susceptible and infected prey, J. Appl. Math. Comput. 284 (2016) 268–285. [5] K. Chakraborty, K. Das, S. Haldar, T.K. Kar, A mathematical study of an eco-epidemiological system on disease persistence and extinction perspective, J. Appl. Math. Comput. 254 (2015) 99–112. [6] B. Sahoo, S. Poria, Diseased prey predator model with general Holling type interactions, J. Appl. Math. Comput. 226 (2014) 83–100. [7] H.W. Hethcote, W. Wang, L. Han, Z. Ma, A predator-prey model with infected prey, J. Theor. Biol. 66 (2004) 259–268. [8] M. Haque, J. Zhen, E. Venturino, An ecoepidemiological predator-prey model with standard disease incidence, J. Math. Methods Appl. Sci. 32 (7) (2009) 875–898. [9] Y. Xiao, L. Chen, Modelling and analysis of a predator-prey model with disease in the prey, J. Math. Biosci. 171 (2001) 59–82. [10] B. Sahoo, S. Poria, Effects of additional food on an ecoepidemic model with time delay on infection, J. Appl. Math. Comput. 245 (2014) 17–35. [11] B. Sahoo, Role of additional food in eco-epidemiological system with diseased in the prey, J. Appl. Math. Comput. 259 (2015) 61–79. [12] P.J. Pal, M. Haque, P.K. Mandal, Dynamics of a predator-prey model with disease in the predator, J. Math. Methods Appl. Sci. 37 (2014) 2429–2450. [13] L. Han, Z. Ma, H.W. Hethcote, Four predator prey models with infectious diseases, J. Math. Comput. Model. 34 (7–8) (2001) 849–858. [14] E. Venturino, Epidemics in predator-prey models: disease in the predators, J. Math. Med. Biol. 19 (3) (2002) 185–205.

Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034

JID: AMC

ARTICLE IN PRESS W. Mbava et al. / Applied Mathematics and Computation 000 (2016) 1–23

[m3Gsc;October 28, 2016;20:34] 23

[15] M.G.L. Mills, L.S. Broomhall, J.T.d. Toit, Cheetah Acinonyx jubatus feeding ecology in the Kruger National Park and a comparison across African Savanna habitats: is the cheetah only a successful hunter on open grassland plains? J. Wildl. Biol. 10 (2004) 177–186. [16] F.G.T. Radloff, J.T.d. Toit, Large predators and their prey in a Southern African Savanna: a predator’s size determines its prey size range, J. Anim. Ecol. 73 (2004) 410–423. [17] C.J. Tambling, S.D. Laurence, S.E. Bellan, E.Z. Cameron, J.T.d. Toit, W.M. Getz, Estimating carnivorian diets using a combination of carcass observations and scats from GPS clusters, J. Zool. 286 (2) (1987) 102–109. [18] Red Data Book of the Mammals of South Africa: A Conservation Assessment, IUCN (World Conservation Union) (2004). [19] M. Maas, D.F. Keet, V.P.M. Rutten, J.A.P. Heesterbeek, M. Nielen, Assessing the impact of feline immunodeficiency virus and bovine tuberculosis co-infection in African lions, Proc. R. Soc. B 279 (2012) 4206–4214. [20] A.R. Renwick, P.C.L. White, R.G. Bengis, Bovine tuberculosis in Southern African wildlife: a multi-species host-pathogen system, J. Epidemiol. Infect. 135 (4) (2007) 529–540. [21] D.F. Keet, P. Hunter, R.G. Bengis, A. Bastos, G.R. Thomson, The 1992 foot-and-mouth disease epizootic in the Kruger National Park, S. Afr. Vet. Assoc. 67 (1996) 83–87. [22] D.F. Keet, A. Michael, D.G.A. Meltzer, Tuberculosis in free-ranging lions (Panthera leo) in the Kruger National Park, in: Proceedings of the South African Veterinary Association Biennial Congress, Durban, Kwazulu-Natal, 20 0 0, pp. 232–241. 20–22 September. [23] Lion (Panthera leo) Bovine Tuberculosis Disease Risk Assessment, Workshop Report, South African National Parks, Endangered Wildlife Trust and Conservation Breeding Specialist Group, 2009. [24] G.T. Skalski, J.F. Gilliam, Functional responses with predator interference: viable alternatives to the Holling type II model, J. Ecol. 82 (2001) 3083–3092. [25] O. Osuna, G. Villaseñor, On the Dulac functions, J. Qual. Theory Dyn. Syst. 10 (1) (2011) 43–49. [26] W.E. Boyce, R.C. DiPrima, Elementary Differential Equations and Boundary Value Problems, ninth ed., Wiley, 2009. [27] A. Hurwitz, et al., On the conditions under which an equation has only roots with negative real parts, in: R. Ballman (Ed.), Report in Selected Papers on Mathematical Trends in Control Theory, Dover, New York, 1964. [28] R.A. Smith, Some applications of Hausdorff dimension inequalities for ordinary differential equations, Proc. R. Soc. Edinb. Sect. A 104 (3–4) (1986) 235–259. [29] M.Y. Li, J. Muldowney, On Bendixson’s criterion, J. Differ. Equ. 106 (1993) 27–39. [30] B. Buonomo, D. Lacitignola, On the use of the geometric approach to global stability for three dimensional ODE systems: a bilinear case, J. Math. Anal. Appl. 348 (1) (2008) 255–266. [31] M.Y. Li, H.L. Smith, L. Wang, Global dynamics of an SEIR epidemic model with vertical transmission, SIAM J. Appl. Math. 62 (1) (2001) 58–69. [32] J.P. Tian, J. Wang, Global stability for cholera epidemic models, J. Math. Biosci. 232 (2011) 31–41. [33] S.B. Hsu, A survey of constructing Lyapunov functions for mathematical models in population biology, Taiwan. J. Math. 9 (2) (2005) 151–173. [34] B. Dubey, R.K. Upadhyay, Persistence and extinction of one-prey and two-predators system, J. Nonlinear Anal. Model. Control 9 (4) (2004) 307–329. [35] H.R. Thieme, Persistence under relaxed point-dissipativity with application to an endemic model, SIAM. J. Math. Anal. 24 (1993) 407–435. [36] Sanparks, Kruger National Park: Biodiversity Statistics. http://www.sanparks.org/parks/kruger/conservation/scientific/ff/biodiversity_statistics.php, 2015, (accessed 04.12.15). [37] D.F. Keet, N.P. Kriek, A. Michael, Tuberculosis and its geographical distribution in free-ranging lions in the Kruger National Park, in: Proceedings of the Third International Conference on Mycobacterium bovis, Cambridge, UK, 20 0 0. [38] K. Marnewick, S.M. Ferreira, S. Grange, J. Watermeyer, N. Maputla, H.T. Davies-Mostert, Evaluating the status of and African wild dogs Lycaon pictus and cheetahs Acinonyx jubatus through tourist-based photographic surveys in the Kruger National Park, PLOS ONE 9 (1) (2014) 1–8. [39] N. Fairhall, Production parameters of the impala, Aepyceros melampus, S. Afr. J. Anim. Sci. 13 (3) (1983) 176–179. [40] Cheetah Populationand Habitat Viability Assessment, Workshop Report, Endangered Wildlife Trust and Conservation Breeding Specialist Group, 2009. [41] M.J. Kelly, S.M. Durant, Viability of the serengeti cheetah population, J. Conserv. Biol. 14 (3) (20 0 0) 786–797. [42] L.S. Broomhall, Cheetah Acinonyx jubatus ecology in the Kruger National Park: A Comparison with other Studies across the Grassland-Woodland Gradient in African Savannas, (M.Sc. thesis), vol. 115, University of Pretoria, South Africa. Unpublished. [43] J.D. Skinner, R.H.N. Smithers, in: The Mammals of Southern African Subregion, University of Pretoria, Pretoria, 1990. [44] P.J. Funston, M.G.L. Mills, The influence of lion predation on the population dynamics of common large ungulates in the Kruger National Park, J. Wildl. Res. 36 (1) (2005) 9–22. [45] J.D.C. Linnell, O. Strand, Interference interactions, co-existence and conservation of mammalian carnivores, J. Divers. Distrib. 6 (20 0 0) 169–176. [46] M.E. Soule, D.T. Bolger, A.C. Alberts, J. Wright, M. Sorice, S. Hill, Reconstructed dynamics of rapid extinctions of chaparral-requiring birds in urban habitat islands, J. Conserv. Biol. 2 (1988) 75–92. [47] P.D. Vickery, M.L. Hunter, J.V. Wells, Evidence of incidental nest predation and its effects on nests of threatened grassland birds, Oikos 63 (1992) 281–288. [48] L. Perko, Differential Equations and Dynamical Systems, third ed., Springer, 2001.

Please cite this article as: W. Mbava et al., Prey, predator and super-predator model with disease in the super-predator, Applied Mathematics and Computation (2016), http://dx.doi.org/10.1016/j.amc.2016.10.034