PARADIGMA BARU PENDIDIKAN MATEMATIKA DAN APLIKASI ONLINE INTERNET PEMBELAJARAN


 
How to cite: A. L. Firdiansyah, “Dynamics of Infected Predator-Prey System with Nonlinear Incidence 

Rate and Prey in Refuge”, J. Mat. Mantik, vol. 6, no. 2, pp. 123-134, October 2020. 

 

Dynamics of Infected Predator-Prey System with Nonlinear 

Incidence Rate and Prey in Refuge  
 

 

Adin Lazuardy Firdiansyah 

STAI Muhammadiyah Probolinggo, adin.lazuardy@gmail.com 

 
doi: https://doi.org/10.15642/mantik.2020.6.2.123-134 

 

 
Abstrak. Model mangsa-pemangsa dengan laju kejadian nonlinear dan perlindungan dimunculkan 

untuk menggambarkan perubahan prilaku pada mangsa yang sehat ketika jumlah mangsa yang 

terinfeksi meningkat, sedangkan pemangsa dapat memakan mangsa dengan mengakses perlindungan 

pada mangsa. Oleh karena itu, analisis dinamik model mangsa-pemangsa dengan penyebaran 

penyakit yang dinotasikan dengan laju kejadian nonlinear dan perlindungan mangsa dibahas pada 

penilitian ini. Hasil analisis ditemukan bahwa ada delapan titik kesetimbangan, dimana semua titik 

kesetimbangan tersebut stabil asimtotik secara lokal. Selanjutnya, perlindungan mangsa juga 

dipertimbangkan pada penelitian ini. Perlindungan mangsa memiliki pengaruh penting pada model. 

Hal ini ditunjukkan dengan adanya perlindungan mangsa yang dapat mencegah kepunahan pada 

populasi mangsa. Pada bagian akhir, simulasi numerik dilakukan untuk mengilustrasikan hasil analisis 

yang diperoleh. Untuk penelitian berikutnya, model mangsa-pemangsa dapat diinvestigasi efek 

pemanenanya untuk kedua populasi. 

 
Kata kunci: Model mangsa-pemangsa; Laju kejadian nonlinear; Perlindungan; Kestabilan 

 
Abstract. A predator-prey system with nonlinear incidence rate and refuging in prey is proposed to 

describe behavior change of certain infected diseases on healthy prey when the number of infected 

prey is getting large, while predator can predate prey by accessing refuging in prey. Therefore, this 

paper discusses the dynamics behavior predator-prey model with the spread of infected disease that 

is denoted by nonlinear incidence rate and adding prey refuge. We find the existence of eight non-

negative equilibrium in the model, which their local stability has been determined. Furthermore, we 

also observe the prey refuge properties in the model. We find that prey refuge can prevent extinction 

in prey populations. In the end, some numerical solutions are carried out to illustrate our analytic 

results. For future work, we can investigate the harvesting effect in both populations, which is disease 

control in the predator-prey model with the spread of infected disease.   

 
Keywords: Predator-prey system; Nonlinear incidence rate; Refuge; Stability 

 

  

Jurnal Matematika MANTIK 

Vol. 6, No. 2, October 2020, pp. 123-134 

 

ISSN: 2527-3159 (print) 2527-3167 (online) 



Jurnal Matematika MANTIK 

Volume 6, No. 2, October 2020, pp. 123-134 

124 

1. Introduction 
 

After Lotka and Volterra proposed a mathematical model that could represent the 

interactions of predators and prey, many researchers were interested to study the model. 

Similarly, Kermack and McKendrick [1] also proposed a SIR (Susceptible-Infectious-

Recovered) epidemic model that represented the infected disease transmission. After that, 

Anderson and May [2] for the first time proposed the eco-epidemiology model and 

investigated invasion, persistence, and spread of infected disease in the model. Since then, 

many researchers have studied the eco-epidemiology model. Generally speaking, the eco-

epidemiology model is divided into three cases. The first case is an only infected disease 

in prey as in [3]–[5]. The second case is an only infected disease in predator as in [6]. 

Meanwhile, the third case is an infected disease in both populations as in [7]–[9].  

In the eco-epidemiology model, there is a basic component that can represent the 

spread of infected disease. Generally, it is denoted by the simple mass incidence rate, 𝛽𝑆𝐼, 
with 𝛽 means the infection rate. This incidence rate shows that the transmission disease 
with rate 𝛽𝑆𝐼, i.e. when the transmission disease increases significantly, then the number 
of infected population increases too [1]. In the model [3]–[8], the authors assume the simple 

mass incidence rate to represent the transmission disease in the model.     

In several cases, the simple mass incidence rate doesn’t produce appropriate results. 

When the amount of infected population increases significantly, the susceptible population 

tend to change their behavior to reduce contact with the infected population. Therefore, 

Capasso and Serio [10] proposed the saturated incidence rate to describe the transmission 

disease. As done in [9], they investigated the local stability where the transmission disease 

in the eco-epidemiology model followed the saturated incidence rate. However, there is a 

nonlinear incidence rate that is suggested by [11]. They introduced nonlinear incidence 

rate, 
𝛽𝐼𝑝

1+𝛼𝐼𝑞
 (𝑝, 𝑞, 𝛼, 𝛽 were positive constants), where 𝛽𝐼𝑝 meant infection rate and 

1

1+𝛼𝐼𝑞
 

meant inhibition rate from the behavior change on healthy population when infected 

population was getting large. Many researcher uses this incidence rate by giving 𝑝, 𝑞, 𝛼 
makes different values as in [12]–[14]. This becomes more rational because it includes the 

effect of tightness from infective individuals [12]. 

Recently, the effect of disease in the eco-epidemiology model has become an 

important topic for many researchers. However, other components that can influence the 

dynamic of species interactions is the Allee effect, habitat complexity, harvesting, and prey 

refuge. For the prey refuge effect, theoretical research provides a conclusion that it can 

influence stabilizing and destabilizing the predator-prey model and can avoid extinction in 

prey [15]–[17]. It can make the predator-prey model to form more realistic.   

Motivated by previous research, we modify a model from [9] by including nonlinear 

incidence rate and adding the effect of prey refuge. At this time, several similar models 

have emerged, but the latest distinctive feature of our model is the inclusion of transmission 

disease in both populations and also the inclusion of prey refuge properties of the prey 

population. Incorporation of prey refuge gives a factor which can be accessed by predator 

populations. Under this adding effect, our model is more realistic and differs from the 

previous eco-epidemiology model. The model is analyzed to determine the local stability 

of its equilibrium points. Moreover, several simulations are given to demonstrate the results 

of our analysis.  
 

2. Methods 
 

The method used for this research is the literature study which studies some previous 

research. It’s used to modify a model from [9]. The steps to solve this research are below: 

a. Reviewing and studying the predator-prey model from previous research. 



Adin Lazuardy Firdiansyah 

Dynamics of Infected Predator-Prey System with Nonlinear Incidence Rate and Prey in Refuge 

125 

b. Modifying the predator-prey model by including nonlinear incidence rate and adding 
the effect of prey refuge.  

c. Analyzing the existence of equilibrium points and local stability in the modified model. 
d. Performing numerical simulation by using the Runge-Kutta 4th order method as a 

numerical method to support our analysis results. 

e. Making a conclusion based on the analysis results.  
 

3. Result and Discussion 
 

3.1 The Formulation Model 
 

Our mathematical model contains two populations, namely the prey population and 

predator population. At time 𝑡, susceptible and infected prey are denoted by 𝑆(𝑡) and 𝐼(𝑡), 
respectively. Meanwhile, susceptible and infected predator are denoted by 𝑌𝑆(𝑡) and 𝑌𝐼(𝑡), 
respectively. The following is several basic assumptions for our mathematical model: 

a. In the absence of predation and disease, the prey grows according to logistics with a 
growth rate 𝑟 (𝑟 > 0) and carrying capacity 𝐾 (𝐾 > 0). 

b. Only one population that can reproduce, namely susceptible prey. 
c. The prey has one infection source like viruses or other sources. Meanwhile, the 

predator can be infected due they eat the infected prey. 

d. The transmission disease follows the nonlinear incidence rate 
𝛽𝑆𝐼

1+𝐼
, where 𝛽𝐼 means 

infection rate and 
1

1+𝐼
 means inhibition rate from behavior change on healthy prey, 

when infected prey is getting large. Meanwhile, the transmission disease in predator 

follows the simple mass incidence rate (𝛾𝑌𝑆𝑌𝐼), where 𝛾 expressed infection rate.   
e. Both of infected prey and predator isn’t recovered and no get immune. 
f. The infected predator can’t predate healthy prey. Both of infected predator and 

susceptible predator can eat infected prey because it is easy to be predated by them. 

But, the ability to catch from an infected predator is lower than susceptible predator. 

g. All species have natural death rates and death rates due to infection.  
h. The functional response for predation of predator is the Lotka-Volterra type.  
i. The refuge protection can be denoted with (1 − 𝑚3)𝑆 for susceptible prey and 

(1 − 𝑚4)𝐼 for infected prey, where 𝑚3, 𝑚4 ∈ (0,1] and healthy prey is more agile 
than infected prey.  

The details of the model structure are shown in the schematic flow diagram as in figure 1. 

From the flow chart in figure 1, the mathematical model is presented as follows: 

𝑑𝑆

𝑑𝑡
= 𝑟𝑆 (1 −

𝑆 + 𝐼

𝐾
) −

𝛽𝑆𝐼

1 + 𝐼
− 𝑝1(1 − 𝑚3)𝑆𝑌𝑆 − 𝑎1𝑆, 

𝑑𝐼

𝑑𝑡
=

𝛽𝑆𝐼

1 + 𝐼
− 𝑝2(1 − 𝑚4)𝐼𝑌𝑆 − 𝑝3(1 − 𝑚4)𝐼𝑌𝐼 − 𝑎2𝐼, 

𝑑𝑌𝑆
𝑑𝑡

= 𝑎3(1 − 𝑚3)𝑆𝑌𝑆 + 𝑎4(1 − 𝑚4)𝐼𝑌𝑆 − 𝛾𝑌𝑆𝑌𝐼 − 𝑑3𝑌𝑆, 

𝑑𝑌𝐼
𝑑𝑡

= 𝑎5(1 − 𝑚4)𝐼𝑌𝐼 + 𝛾𝑌𝑆𝑌𝐼 − 𝑎6𝑌𝐼. 

(1) 

with the initial condition as 𝑆(0) > 0, 𝐼(0) > 0, 𝑌𝑆(0) > 0, 𝑌𝐼(0) > 0. All parameters with 
their biological meaning are given as follows: 𝑎1 = 𝑚1 + 𝑑1, 𝑎2 = 𝑚2 + 𝑑2 + 𝑐, 𝑎3 =
𝑝1𝑞1, 𝑎4 = 𝑝2𝑞2, 𝑎5 = 𝑝3𝑞3, and 𝑎6 = 𝑑4 + 𝑑5 where 𝑚1 and 𝑚2 are migration rates for 
susceptible prey and infected prey, respectively. 𝑝1(𝑝2) and 𝑝3 are predation rate from 
healthy predator to healthy prey (to infected prey) and predation rate from infected predator 

to infected prey, respectively. 𝑞1(𝑞2) and 𝑞3 are conversion rate from healthy prey into 
healthy predator (infected predator) and conversion rate from infected prey into an infected 



Jurnal Matematika MANTIK 

Volume 6, No. 2, October 2020, pp. 123-134 

126 

predator. 𝑑1(𝑑2) and 𝑑3(𝑑4) are natural death from susceptible prey (infected prey) and 
natural death from a susceptible predator (infected predator), respectively. 𝑐(𝑑5) is the 
death rate from infected prey (from an infected predator) due to infection from disease. We 

assume that all parameters are positive values. 

 

Figure 1. Schematic diagram of the model 

 

3.2 Equilibrium Point of Mathematical Model  
 

We set the right-hand sides equal to zero. Thus, we get possible equilibriums:   

a. The trivial equilibrium point is 𝐸0(0,0,0,0). This point always exists.  

b. The axial equilibrium point is 𝐸1(𝑆
(1), 0,0,0), where 𝑆(1) =

𝐾(𝑟−𝑎1)

𝑟
. It will exist when 

𝑟 > 𝑎1.  

c. 𝐸2(𝑆
(2), 𝐼(2), 0,0) is the predator-free equilibrium point, where  

𝑆(2) =
(1 + 𝐼(2))(𝑟𝐾 − 𝑟𝐼(2) − 𝐾𝑎1) − 𝐾𝛽𝐼

(2)

𝑟(1 + 𝐼(2))
 

and 𝐼(2) is the positive root of quadratic equations 𝐴0(𝐼
(2))

2
+ 𝐴1𝐼

(2) + 𝐴2 = 0, 

where 𝐴0 = 𝑟(𝑎2 + 𝛽), 𝐴1 = 2𝑟𝑎2 − 𝛽𝑟𝐾 + 𝛽𝑟 + 𝛽𝑎1𝐾 + 𝛽
2𝐾, and 𝐴2 = 𝑟𝑎2 −

𝛽𝑟𝐾 + 𝛽𝑎1𝐾. The equilibrium 𝐸2 will exist if 𝑟 ≥
𝛽𝑎1𝐾

𝛽𝐾−𝑎2
. 

d. 𝐸3 (𝑆
(3), 0, 𝑌𝑆

(3)
, 0) is the disease-free equilibrium point, where  

𝑆(3) =
𝑑3

𝑎3(1 − 𝑚3)
 

and 

𝑌𝑆
(3)

=
𝐾(𝑟 − 𝑎1) − 𝑟𝑆

(3)

𝐾𝑝1(1 − 𝑚3)
. 

The equilibrium point 𝐸3 will exist if 𝑟 >
𝑎1𝑎3𝐾(1−𝑚3)

𝑎3𝐾(1−𝑚3)−𝑑3
.  

e. 𝐸4 (𝑆
(4),

𝑎6

𝑎5(1−𝑚4)
, 0, 𝑌𝐼

(4)
) is the healthy predator-free equilibrium point, where 

𝑆 𝐼 

𝑌𝑆 
𝑌𝐼 

𝛽𝑆𝐼

1 + 𝐼
 

𝑟𝑆 (1 −
𝑆 + 𝐼

𝐾
) 

𝑎1𝑆 

𝑝1𝑆𝑌𝑆 

𝑑3𝑌𝑆 

𝑚3𝑆𝑌𝑆  

𝑝2𝐼𝑌𝑆 

𝑚4𝐼𝑌𝑆 

𝑝3𝐼𝑌𝐼 

𝛾𝑌𝑆𝑌𝐼  

𝑎2𝐼 

𝑚4𝐼𝑌𝐼 

𝑎6𝑌𝐼  



Adin Lazuardy Firdiansyah 

Dynamics of Infected Predator-Prey System with Nonlinear Incidence Rate and Prey in Refuge 

127 

𝑆(4) =
(𝑎5(1 − 𝑚4) + 𝑎6)(𝑝3(1 − 𝑚4)𝑌𝐼

(4)
+ 𝑎2)

𝛽𝑎5(1 − 𝑚4)
 

and 

𝑌𝐼
(4)

=
𝐶1 − 𝐶2

𝑟𝑝3(1 − 𝑚4)(𝑎5(1 − 𝑚4) + 𝑎6)
2
 

with 𝐶1 = 𝐾𝛽𝑎5(1 − 𝑚4)(𝑟 − 𝑎1)(𝑎5(1 − 𝑚4) + 𝑎6) and 𝐶2 = 𝑟𝑎2(𝑎5(1 − 𝑚4) +
𝑎6)

2 + 𝐾𝛽2𝑎5𝑎6(1 − 𝑚4) + 𝑟𝛽𝑎6(𝑎5(1 − 𝑚4) + 𝑎6). This equilibrium exists when 
𝐶1 > 𝐶2 and 𝑟 > 𝑎1.  

f. 𝐸5 (𝑆
(5), 0,

𝑎6

𝛾
, 𝑌𝐼

(5)
) is the infected prey-free equilibrium point, where  

𝑆(5) =
𝐾𝛾(𝑟 − 𝑎1) − 𝐾𝑝1𝑎6(1 − 𝑚3)

𝑟𝛾
 

and 

𝑌𝐼
(5)

=
(1 − 𝑚3)𝑎3𝑆

(5) − 𝑑3
𝛾

. 

It will exist if 𝑟 >
𝐾𝑎3(1−𝑚3)(𝛾𝑎1+𝑝1𝑎6(1−𝑚3))

𝛾(𝐾𝑎3(1−𝑚3)−𝑑3)
.   

g. 𝐸6 (
𝑑3−𝑎4(1−𝑚4)𝐼

(6)

𝑎3(1−𝑚3)
, 𝐼(6), 𝑌𝑆

(6)
, 0) is the infected predator-free equilibrium, where  

𝑌𝑆
(6)

=
𝛽𝑑3 − 𝛽𝑎4(1 − 𝑚4)𝐼

(6) − 𝑎2𝑎3(1 − 𝑚3)(1 + 𝐼
(6))

𝑝2𝑎3(1 − 𝑚3)(1 − 𝑚4)(1 + 𝐼
(6))

 

and 𝐼(6) is the positive root of quadratic equations 𝐷0(𝐼
(6))

2
+ 𝐷1𝐼

(6) + 𝐷2 = 0 with 

𝐷0 = 𝑟𝑝2(1 − 𝑚4)(𝑎3(1 − 𝑚3) − 𝑎4(1 − 𝑚4)), 
𝐷1 = 𝐾𝑎3𝑝2(1 − 𝑚3)(1 − 𝑚4)(𝛽 + 𝑎1 − 𝑟) 
 +𝑟𝑝2(1 − 𝑚4)(𝑑3 + 𝑎3(1 − 𝑚3) − 𝑎4(1 − 𝑚4)) 

 −𝐾𝑝1(1 − 𝑚3)(𝑎2𝑎3(1 − 𝑚3) + 𝛽𝑎4(1 − 𝑚4)), 
𝐷2 = 𝑟𝑝2(1 − 𝑚4)(𝑑3 − 𝐾𝑎3(1 − 𝑚3)) 
 +𝐾𝑎3(1 − 𝑚3)(𝑎1𝑝2(1 − 𝑚4) − 𝑎2𝑝1(1 − 𝑚3)) + 𝐾𝛽𝑝1𝑑3(1 − 𝑚3). 

This equilibrium will exist if 𝐼(6) <
𝛽𝑑3−𝑎2𝑎3(1−𝑚3)

𝛽𝑎4(1−𝑚4)+𝑎2𝑎3(1−𝑚3)
.   

 

h. 𝐸7 (𝑆
(7), 𝐼(7),

𝑎6−𝑎5(1−𝑚4)𝐼
(7)

𝛾
,
𝑎3(1−𝑚3)𝑆

(7)+𝑎4(1−𝑚4)𝐼
(7)−𝑑3

𝛾
) is the interior equilibrium 

point, where 𝑆(7) =
(1+𝐼(7))(𝐾𝑟𝛾−𝑎1𝛾𝐾−𝑟𝛾𝐼

(7)−𝐾𝑝1(1−𝑚3)(𝑎6−𝑎5(1−𝑚4)𝐼
(7)))−𝐾𝛾𝛽𝐼(7)

𝛾𝑟(1+𝐼(7))
 and  

𝐼(7) is the positive root of the cubic equation 𝑄0(𝐼
(7))

3
+ 𝑄1(𝐼

(7))
2
+ 𝑄2𝐼

(7) + 𝑄3 =

0, with 

𝑄0 = 𝐾𝑎3𝑎5𝑝1𝑝3(1 − 𝑚3)
2(1 − 𝑚4)

2 − 𝛾𝑟𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4) 



Jurnal Matematika MANTIK 

Volume 6, No. 2, October 2020, pp. 123-134 

128 

 +𝛾𝑟𝑎4𝑝3(1 − 𝑚4)
2 − 𝛾𝑟𝑎5𝑝2(1 − 𝑚4)

2, 
𝑄1 = 2𝛾𝑟𝑎4𝑝3(1 − 𝑚4)

2 − 2𝛾𝑟𝑎5𝑝2(1 − 𝑚4)
2 + 𝛾𝑟𝑎6𝑝2(1 − 𝑚4) + 𝛾

2𝑟𝑎2 
 −𝛾𝑟𝑑3𝑝3(1 − 𝑚4) + 2𝐾𝑎3𝑎5𝑝1𝑝3(1 − 𝑚3)

2(1 − 𝑚4)
2 + 𝛽𝛾2𝑟 

 −𝐾𝑎3𝑎6𝑝1𝑝3(1 − 𝑚3)
2(1 − 𝑚4) − 𝐾𝛾𝑎1𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4) 

 +𝐾𝛾𝑟𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4) − 𝐾𝛽𝛾𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4) 
 −𝐾𝛽𝛾𝑎5𝑝1(1 − 𝑚3)(1 − 𝑚4) − 2𝛾𝑟𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4), 
𝑄2 = 𝛾𝑟𝑎4𝑝3(1 − 𝑚4)

2 − 𝛾𝑟𝑎5𝑝2(1 − 𝑚4)
2 + 2𝛾𝑟𝑎6𝑝2(1 − 𝑚4) + 2𝛾

2𝑟𝑎2 
 −2𝛾𝑟𝑑3𝑝3(1 − 𝑚4) + 𝐾𝑎3𝑎5𝑝1𝑝3(1 − 𝑚3)

2(1 − 𝑚4)
2 + 𝛽𝛾2𝑟 

 −2𝐾𝑎3𝑎6𝑝1𝑝3(1 − 𝑚3)
2(1 − 𝑚4) − 2𝐾𝛾𝑎1𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4) 

 +2𝐾𝛾𝑟𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4) − 𝐾𝛽𝛾𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4) + 𝐾𝛽
2𝛾2 

 −𝐾𝛽𝛾𝑎5𝑝1(1 − 𝑚3)(1 − 𝑚4) + 𝐾𝛽𝛾𝑎6𝑝1(1 − 𝑚3) − 𝐾𝛽𝛾
2𝑟 + 𝐾𝛽𝛾2𝑎1 

 −𝛾𝑟𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4), 
𝑄3 = 𝐾𝛾𝑟𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4) − 𝐾𝛾𝑎1𝑎3𝑝3(1 − 𝑚3)(1 − 𝑚4) 
 −𝐾𝑎3𝑎6𝑝1𝑝3(1 − 𝑚3)

2(1 − 𝑚4) + 𝐾𝛽𝛾𝑎6𝑝1(1 − 𝑚3) 

 
+𝛾𝑟𝑎6𝑝2(1 − 𝑚4) − 𝛾𝑟𝑑3𝑝3(1 − 𝑚4) + 𝐾𝛽𝛾

2𝑎1 − 𝐾𝛽𝛾
2𝑟 + 𝛾2𝑟𝑎2. 

The existing form of the positive root in the cubic equation can be determined by using 

Cardan’s method as in [18]. 
 

3.3 Stability Analysis 
 

To check the local stability of the equilibrium point, we have to determine the 

eigenvalues of the Jacobian matrix. Here, the Jacobian matrix at the equilibrium point 

𝐸∗(𝑆∗, 𝐼∗, 𝑌𝑆
∗, 𝑌𝐼

∗) is as follows: 

𝐽∗ =

[
 
 
 
 

𝑢11 𝑢12 −𝑝1(1 − 𝑚3)𝑆
∗ 0

𝛽𝐼∗

1 + 𝐼∗
𝑢22 −𝑝2(1 − 𝑚4)𝐼

∗ −𝑝3(1 − 𝑚4)𝐼
∗

𝑎3(1 − 𝑚3)𝑌𝑆
∗ 𝑎4(1 − 𝑚4)𝑌𝑆

∗ 𝑢33 −𝛾𝑌𝑆
∗

0 𝑎5(1 − 𝑚4)𝑌𝐼
∗ 𝛾𝑌𝐼

∗ 𝑢44 ]
 
 
 
 

, 

where 

𝑢11 = = 𝑟 (1 −
𝑆∗ + 𝐼∗

𝐾
) −

𝑟𝑆∗

𝐾
−

𝛽𝐼∗

1 + 𝐼∗
− 𝑝1(1 − 𝑚3)𝑌𝑆

∗ − 𝑎1, 

𝑢12 = −
𝑟𝑆∗

𝐾
−

𝛽𝑆∗

1 + 𝐼∗
+

𝛽𝑆∗𝐼∗

(1 + 𝐼∗)2
, 

𝑢22 = 
𝛽𝑆∗

1 + 𝐼∗
−

𝛽𝑆∗𝐼∗

(1 + 𝐼∗)2
− 𝑝2(1 − 𝑚4)𝑌𝑆

∗ − 𝑝3(1 − 𝑚4)𝑌𝐼
∗ − 𝑎2, 

𝑢33 = 𝑎3(1 − 𝑚3)𝑆
∗ + 𝑎4(1 − 𝑚4)𝐼

∗ − 𝛾𝑌𝐼
∗ − 𝑑3, 

𝑢44 = 𝑎5(1 − 𝑚4)𝐼
∗ + 𝛾𝑌𝐼

∗ − 𝑎6. 

The characteristic equation of 𝐸0 is  

(𝑟 − 𝑎1 − 𝜆)(−𝑎2 − 𝜆)(−𝑑3 − 𝜆)(−𝑎6 − 𝜆) = 0, 

which has eigenvalues 𝜆1 = 𝑟 − 𝑎1; 𝜆2 = −𝑎2; 𝜆3 = −𝑑3; 𝜆4 = −𝑎6. Thus, 𝐸0 is locally 
asymptotically stable if only if 𝑟 < 𝑎1. On point 𝐸1, the characteristic equation of 𝐸1 is  

(𝑎1 − 𝑟 − 𝜆)(𝑔1 − 𝜆)(𝑔2 − 𝜆)(−𝑎6 − 𝜆) = 0, 

with 



Adin Lazuardy Firdiansyah 

Dynamics of Infected Predator-Prey System with Nonlinear Incidence Rate and Prey in Refuge 

129 

𝑔1 = 
𝛽𝐾(𝑟 − 𝑎1)

𝑟
− 𝑎2, 

𝑔2 = 
𝑎3𝐾(1 − 𝑚3)(𝑟 − 𝑎1)

𝑟
− 𝑑3. 

This point has eigenvalues 𝜆1 = 𝑎1 − 𝑟; 𝜆2 = 𝑔1; 𝜆3 = 𝑔2; 𝜆4 = −𝑎6. Thus, 𝐸1 can be said 

the locally asymptotically stable if only if 𝑎1 < 𝑟 < 𝑚𝑖𝑛 {
𝛽𝐾𝑎1

𝛽𝐾−𝑎2
,

𝐾𝑎1𝑎3(1−𝑚3)

𝐾𝑎3(1−𝑚3)−𝑑3
}, where 

𝛽𝐾 > 𝑎2 and 𝐾𝑎3(1 − 𝑚3) > 𝑑3. On point 𝐸2, the characteristic equation of 𝐸2 is  

(𝜆2 − (ℎ1 + ℎ2)𝜆 + (ℎ1ℎ2 − ℎ3ℎ2))(ℎ4 − 𝜆)(ℎ5 − 𝜆) = 0, 

with 

ℎ1 = 𝑟 (1 −
𝑆(2) + 𝐼(2)

𝐾
) −

𝑟𝑆(2)

𝐾
−

𝛽𝐼(2)

1 + 𝐼(2)
− 𝑎1, 

 

ℎ2 = 
𝛽𝑆(2)

1 + 𝐼(2)
−

𝛽𝑆(2)𝐼(2)

(1 + 𝐼(2))2
− 𝑎2, 

 

ℎ3 = 
𝛽𝐼(2)

1 + 𝐼(2)
, 

 

ℎ4 = 𝑎3(1 − 𝑚3)𝑆
(2) + 𝑎4(1 − 𝑚4)𝐼

(2) − 𝑑3,  

ℎ5 = 𝑎5(1 − 𝑚4)𝐼
(2) − 𝑎6.  

Eigenvalues on this point 𝐸2 is 𝜆1,2 =
(ℎ1+ℎ2)±√(ℎ1+ℎ2)

2−4(ℎ1ℎ2−ℎ3ℎ2)

2
; 𝜆3 = ℎ4; 𝜆4 = ℎ5. 

Thus, the point 𝐸2 can be said the locally asymptotically stable if only if ℎ1 + ℎ2 < 0; 
ℎ1ℎ2 > ℎ3ℎ2; 𝑎5(1 − 𝑚4)𝐼

(2) < 𝑎6; 𝑎3(1 − 𝑚3)𝑆
(2) + 𝑎4(1 − 𝑚4)𝐼

(2) < 𝑑3. On point 
𝐸3,  we obtain the characteristic equation of 𝐸3 is  

(𝜆2 − (𝑛1 + 𝑛2)𝜆 + (𝑛1𝑛2 − 𝑛3𝑛4))(𝑛5 − 𝜆)(𝑛6 − 𝜆) = 0, 

with 

𝑛1 = 𝑟 (1 −
𝑆(3)

𝐾
) −

𝑟𝑆(3)

𝐾
− 𝑝1(1 − 𝑚3)𝑌𝑆

(3)
− 𝑎1, 

 

𝑛2 = 𝑎3(1 − 𝑚3)𝑆
(3) − 𝑑3,  

𝑛3 = −𝑝1(1 − 𝑚3)𝑆
(3),  

𝑛4 = 𝑎3(1 − 𝑚3)𝑌𝑆
(3)

,  

𝑛5 = 𝛽𝑆(3) − 𝑝2(1 − 𝑚4)𝑌𝑆
(3)

− 𝑎2, 
 

𝑛6 = 𝛾𝑌𝑆
(3)

− 𝑎6, 
 

which has eigenvalues 𝜆1,2 =
(𝑛1+𝑛2)±√(𝑛1+𝑛2)

2−4(𝑛1𝑛2−𝑛3𝑛4)

2
; 𝜆3 = 𝑛5; 𝜆4 = 𝑛6. So, we 

can say that the point 𝐸3 is locally asymptotically stable if only if 𝑛1 + 𝑛2 < 0; 

𝑛1𝑛2 > 𝑛3𝑛4; 𝛽𝑆
(3) − 𝑝2(1 − 𝑚4)𝑌𝑆

(3)
< 𝑎2 and 𝛾𝑌𝑆

(3)
< 𝑎6. On point 𝐸4, the 

characteristic equation of 𝐸4 is  

(𝑧1 − 𝜆)(𝜆
3 + 𝜑1𝜆

2 + 𝜑2𝜆 + 𝜑3) = 0. 

One eigenvalue is 𝜆1 = 𝑧1, where 𝑧1 = 𝑎3(1 − 𝑚3)𝑆
(4) + 𝑎4(1 − 𝑚4)𝐼

(4) − 𝛾𝑌𝐼
(4)

− 𝑑3. 
Other eigenvalues are obtained from the roots of cubic equations, which are obtained from 

the matrix 𝐽(𝐸4).  



Jurnal Matematika MANTIK 

Volume 6, No. 2, October 2020, pp. 123-134 

130 

𝐽(𝐸4) =

(

 
 

𝑧2 𝑧3 0

𝛽𝐼(4)

1 + 𝐼(4)
𝑧4 −𝑝3(1 − 𝑚4)𝐼

(4)

0 𝑎3(1 − 𝑚4)𝑌𝐼
(4)

𝑎5(1 − 𝑚4)𝐼
(4) − 𝑎6)

 
 

, 

with  

𝑧2 = 𝑟 (1 −
𝑆(4) + 𝐼(4)

𝐾
) −

𝑟𝑆(4)

𝐾
−

𝛽𝐼(4)

1 + 𝐼(4)
− 𝑎1, 

𝑧3 = −
𝑟𝑆(4)

𝐾
−

𝛽𝑆(4)

1 + 𝐼(4)
+

𝛽𝑆(4)𝐼(4)

(1 + 𝐼(4))2
, 

𝑧4 = 
𝛽𝑆(4)

1 + 𝐼(4)
−

𝛽𝑆(4)𝐼(4)

(1 + 𝐼(4))2
− 𝑝3(1 − 𝑚4)𝑌𝐼

(4)
− 𝑎2. 

Thus, the roots of the cubic equation are 𝜑1 = −𝑡𝑟 (𝐽(𝐸4)) ; 𝜑2 = 𝑀11 + 𝑀22 + 𝑀33; and 

𝜑3 = − 𝑑𝑒𝑡 (𝐽(𝐸4)), where 𝑀𝑖𝑗 is minor of entry (𝐽(𝐸4))
𝑖𝑗
. By using Routh-Hurwitz 

criteria, the point 𝐸4 is the locally asymptotically stable which can be provided with the 

following condition 𝜑1 > 0; 𝜑3 > 0; 𝜑1𝜑2 > 𝜑3; 𝑎3(1 − 𝑚3)𝑆
(4) + 𝑎4(1 − 𝑚4)𝐼

(4) <

𝛾𝑌𝐼
(4)

+ 𝑑3. On point 𝐸5, the characteristic equation of 𝐸5 is  

(𝑒1 − 𝜆)(𝜆
3 + 𝛿1𝜆

2 + 𝛿2𝜆 + 𝛿3) = 0. 

One eigenvalue is 𝜆1 = 𝑒1, where 𝑒1 = 𝛽𝑆
(5) − 𝑝2(1 − 𝑚4)𝑌𝑆

(5)
− 𝑝3(1 − 𝑚4)𝑌𝐼

(5)
− 𝑎2. 

Meanwhile, three eigenvalues are the roots of cubic equations 𝜆3 + 𝛿1𝜆
2 + 𝛿2𝜆 + 𝛿3 = 0, 

which are given by matrix 𝐽(𝐸5).  

𝐽(𝐸5) = (

𝑒2 −𝑝1(1 − 𝑚3)𝑆
(5) 0

𝑎3(1 − 𝑚3)𝑌𝑆
(5)

𝑒3 −𝛾𝑌𝑆
(5)

0 𝛾𝑌𝐼
(5)

𝑒4

), 

where 

𝑒2 = 𝑟 (1 −
𝑆(5)

𝐾
) −

𝑟𝑆(5)

𝐾
− 𝑝1(1 − 𝑚3)𝑌𝑆

(5)
− 𝑎1, 

𝑒3 = 𝑎3(1 − 𝑚3)𝑆
(5) − 𝛾𝑌𝐼

(5)
− 𝑑3, 

𝑒4 = 𝛾𝑌𝑆
(5)

− 𝑎6. 

So, we obtain the roots of the cubic equation, i.e. 𝛿1 = −𝑡𝑟 (𝐽(𝐸5)) ; 𝛿2 = 𝑀11 + 𝑀22 +

𝑀33; and 𝛿3 = − 𝑑𝑒𝑡 (𝐽(𝐸5)) with 𝑀𝑖𝑗 is minor of entry (𝐽(𝐸5))
𝑖𝑗
. By using Routh-

Hurwitz criteria, the point 𝐸5 can be said locally asymptotically stable if  it satisfies the 

following condition, 𝛽𝑆(5) < 𝑝2(1 − 𝑚4)𝑌𝑆
(5)

+ 𝑝3(1 − 𝑚4)𝑌𝐼
(5)

+ 𝑎2; 𝛿1 > 0; 𝛿3 >
0; 𝛿1𝛿2 > 𝛿3. On point 𝐸6, the characteristic equation of 𝐸6 is  

(𝑙1 − 𝜆)(𝜆
3 + 𝜎1𝜆

2 + 𝜎2𝜆 + 𝜎3) = 0. 

One eigenvalue is 𝜆1 = 𝑙1, with 𝑙1 = 𝑎5(1 − 𝑚4)𝐼
(6) + 𝛾𝑌𝑆

(6)
− 𝑎6. Three eigenvalues are 

the roots of cubic equations, i.e. 𝜆3 + 𝜎1𝜆
2 + 𝜎2𝜆 + 𝜎3 = 0, which are given by matrix 



Adin Lazuardy Firdiansyah 

Dynamics of Infected Predator-Prey System with Nonlinear Incidence Rate and Prey in Refuge 

131 

𝐽(𝐸6). 

𝐽(𝐸6) =

(

 
 

𝑙2 𝑙3 −𝑝1(1 − 𝑚3)𝑆
(6)

𝛽𝐼(6)

1 + 𝐼(6)
𝑙4 −𝑝2(1 − 𝑚4)𝐼

(6)

𝑎3(1 − 𝑚3)𝑌𝑆
(6)

𝑎4(1 − 𝑚4)𝑌𝑆
(6)

 𝑙5 )

 
 

, 

where 

𝑙2 = 𝑟 (1 −
𝑆(6) + 𝐼(6)

𝐾
) −

𝑟𝑆(6)

𝐾
−

𝛽𝐼(6)

1 + 𝐼(6)
− 𝑝1(1 − 𝑚3)𝑌𝑆

(6)
− 𝑎1, 

𝑙3 = −
𝑟𝑆(6)

𝐾
−

𝛽𝑆(6)

1 + 𝐼(6)
+

𝛽𝑆(6)𝐼(6)

(1 + 𝐼(6))2
, 

𝑙4 = 
𝛽𝑆(6)

1 + 𝐼(6)
−

𝛽𝑆(6)𝐼(6)

(1 + 𝐼(6))2
− 𝑝2(1 − 𝑚4)𝑌𝑆

(6)
− 𝑎2, 

𝑙5 = 𝑎3(1 − 𝑚3)𝑆
(6) + 𝑎4(1 − 𝑚4)𝐼

(6) − 𝑑3. 

Thus, we obtained other corresponding eigenvalues, namely 𝜎1 = −𝑡𝑟 (𝐽(𝐸6)) ; 

𝜎2 = 𝑀11 + 𝑀22 + 𝑀33; 𝜎3 = − 𝑑𝑒𝑡 (𝐽(𝐸6)), with a minor of entry (𝐽(𝐸6))
𝑖𝑗

 that is 

denoted by 𝑀𝑖𝑗. By using Routh-Hurwitz criteria, the point 𝐸6 can be said locally 

asymptotically stable if only if 𝑎5(1 − 𝑚4)𝐼
(6) + 𝛾𝑌𝑆

(6)
< 𝑎6; 𝜎1 > 0; 𝜎3 > 0; 𝜎1𝜎2 > 𝜎3. 

On point 𝐸7, we have the characteristic equation of 𝐸7, which is obtained from the matrix 
𝐽∗(𝐸7). Meanwhile, 𝐽

∗(𝐸7) is the Jacobian matrix at 𝐸7. Eigenvalues of 𝐸7 are the roots of 
the fourth-order equation, i.e. 

𝜆4 + 𝜂1𝜆
3 + 𝜂2𝜆

2 + 𝜂3𝜆 + 𝜂4 = 0, 

with 

𝜂1 = −𝑡𝑟(𝐽
∗(𝐸7)), 

𝜂2 = 𝜏11 + 𝜏22 + 𝜏33 + 11 + 22 + 𝜌11, 
𝜂4 = −(𝑀11 + 𝑀22 + 𝑀33 + 𝑀44), 
𝜂4 = 𝑑𝑒𝑡(𝐽

∗(𝐸7)), 

where 𝜏( ) and 𝜌 are minor of order 2 from the remaining submatrix after 4th row and 4th 
column (after 3rd row and 3rd column) and after 2nd row and 2nd column are removed from 

𝐽∗(𝐸7), respectively. Meanwhile, 𝑀𝑖𝑗 is minor of order 3 for entry (𝐽
∗(𝐸7))𝑖𝑗

. By using 

Routh-Hurwitz criteria, we get the condition for the stability of asymptotical locally at the 

point 𝐸7, namely 𝜂1 > 0; 𝜂1𝜂2 > 𝜂3; 𝜂3(𝜂2𝜂1 − 𝜂3) − 𝜂4𝜂1
2 > 0; 𝜂4 > 0. 

 

3.4 Numerical Simulation 
 

In this section, we present some numerical solutions by using the Runge-Kutta 4th 
order as a numerical method to illustrate our analytical results. We use some parameters 

for the system (1), namely 𝑟 = 0.45, 𝐾 = 10, 𝛾 = 0.15, 𝑝1 = 0.6, 𝑝2 = 0.8, 𝑝3 =
0.45, 𝑎1 = 0.12, 𝑎2 = 0.16, 𝑎3 = 0.25, 𝑎4 = 0.35, 𝑎5 = 0.15, 𝑎6 = 0.16, 𝑑3 =
0.103, 𝛽 = 0.3333.  
 

Case 1: We start with 𝑚3 = 0.157 and 𝑚4 = 0.108. The prey refuge is very small, which 
means that the number of prey outside the refuge is huge. All solutions are convergent to 



Jurnal Matematika MANTIK 

Volume 6, No. 2, October 2020, pp. 123-134 

132 

𝐸3(0.4887,0,0.6090,0), see figure 2(a). This equilibrium is also locally asymptotically 
stable, see figure 2(b). It indicates that infected predator and infected prey are extinct. 

 

Case 2: We consider 𝑚3 = 0.65 and 𝑚4 = 0.3. The prey refuge means that the number of 
prey outside prey refuge is less than the first case. All solutions with different values are 

convergent to equilibrium 𝐸5(2.3556,0,1.0667,0.6874), see figure 3(a). 𝐸5 is also locally 
asymptotically stable, see figure 3(b). It indicates that the infected prey is extinct. 

  
(a) (b) 

Figure 2. Numerical solution for case 1: (a) time graph; (b) the phase portrait 

  
(a) (b) 

Figure 3. Numerical solution for case 2: (a) time graph; (b) the phase portrait 

 

Case 3: When the prey refuge increases to become 𝑚3 = 0.85 and 𝑚4 = 0.5, then it shows 
that amount of prey outside refuge is small. All trajectories with different values converge 

to the point of equilibrium 𝐸7(2.7641,0.4901,0.8216,0.5761), see figure 4(a), and is also 
locally asymptotically stable, see figure 4(b). It indicates that all populations exist. 

  
(a) (b) 

Figure 4. Numerical solution for case 3: (a) time graph; (b) the phase portrait 

 

To investigate the effect of refuge in prey populations, we can observe the dynamic 

behavior of the system (1) by using 𝑚3, 𝑚4 makes different values. Figure 5 shows the 



Adin Lazuardy Firdiansyah 

Dynamics of Infected Predator-Prey System with Nonlinear Incidence Rate and Prey in Refuge 

133 

time graph of prey populations which means that prey refuge can prevent extinction in prey 

populations when 𝑚3, 𝑚4 values are getting large. 
 

 

 
Figure 5. Influence refuge in prey by using 𝑚3, 𝑚4 makes different values 

   

4.   Conclusions 
 

In this paper, we have merged an eco-epidemiology model with transmission disease 

in both populations by using non-linear incidence rate in prey populations and also prey 

refuge proportional in prey populations. Eight equilibrium points that exist under certain 

conditions and also that local stability for those points have been determined. We have 

noticed that prey refuge can avoid extinction in prey populations. For future work, we have 

to investigate the harvesting effect in both populations. It is used as a disease control in an 

eco-epidemiology model. 
 

References 
 

[1] W. O. Kermack and A. G. McKendrick, “A Contribution to the Mathematical Theory 

of Epidemics,” Proc. R. Soc. London, vol. 115, no. 772, p. 22, 1927. 

[2] R. M. Anderson and R. M. May, “The invasion, persistence and spread of infectious 

diseases within animal and plant communities,” Philos. Trans. R. Soc. Lond. B. Biol. 

Sci., vol. 314, no. 1167, pp. 533–570, 1986, doi: 10.1098/rstb.1986.0072. 

[3] C. Maji, D. Kesh, and D. Mukherjee, “Bifurcation and global stability in an eco-

epidemic model with refuge,” Energy, Ecol. Environ., vol. 4, no. 3, pp. 103–115, 

2019, doi: 10.1007/s40974-019-00117-6. 

[4] S. Kant and V. Kumar, “Analysis of an eco-epidemiological model with migrating 

and refuging prey,” Springer Proc. Math. Stat., vol. 143, pp. 17–36, 2015, doi: 

10.1007/978-81-322-2485-3_2. 

[5] A. K. Pal and G. P. Samanta, “Stability analysis of an eco-epidemiological model 

incorporating a prey refuge,” Nonlinear Anal. Model. Control, vol. 15, no. 4, pp. 473–

491, 2010, doi: 10.15388/na.15.4.14319. 

[6] S. A. Wuhaib and Y. A. B. U. Hasan, “Predator-Prey Interactions With Harvesting 

of Predator With Prey in Refuge,” Commun. Math. Biol. Neurosci., vol. 2013, no. 1, 

pp. 1–19, 2013. 

[7] S. P. Bera, A. Maiti, and G. P. Samanta, “A prey-predator model with infection in 

both prey and predator,” Filomat, vol. 29, no. 8, pp. 1753–1767, 2015, doi: 

10.2298/FIL1508753B. 



Jurnal Matematika MANTIK 

Volume 6, No. 2, October 2020, pp. 123-134 

134 

[8] S. Kant and V. Kumar, “Stability analysis of predator–prey system with migrating 

prey and disease infection in both species,” Appl. Math. Model., vol. 42, pp. 509–

539, 2017, doi: 10.1016/j.apm.2016.10.003. 

[9] A. Pusparani, W. M. Kusumawinahyu, and Trisilowati, “Dynamical Analysis of 

Infected Predator-Prey Model with Saturated Incidence Rate,” IOP Conf. Ser. Mater. 

Sci. Eng., vol. 546, no. 5, p. 8, 2019, doi: 10.1088/1757-899X/546/5/052055. 

[10] V. Capasso and G. Serio, “A generalization of the Kermack-McKendrick 

deterministic epidemic model,” Math. Biosci., vol. 42, no. 1–2, pp. 43–61, 1978, doi: 

10.1016/0025-5564(78)90006-8. 

[11] W. min Liu, S. A. Levin, and Y. Iwasa, “Influence of nonlinear incidence rates upon 

the behavior of SIRS epidemiological models,” J. Math. Biol., vol. 23, no. 2, pp. 187–

204, 1986, doi: 10.1007/BF00276956. 

[12] S. Ruan and W. Wang, “Dynamical behavior of an epidemic model with a nonlinear 

incidence rate,” J. Differ. Equ., vol. 188, no. 1, pp. 135–163, 2003, doi: 

10.1016/S0022-0396(02)00089-X. 

[13] R. K. Naji and A. N. Mustafa, “The dynamics of an eco-epidemiological model with 

nonlinear incidence rate,” J. Appl. Math., vol. 2012, p. 24, 2012, doi: 

10.1155/2012/852631. 

[14] A. P. Maiti, C. Jana, and D. K. Maiti, “A delayed eco-epidemiological model with 

nonlinear incidence rate and Crowley–Martin functional response for infected prey 

and predator,” Nonlinear Dyn., vol. 98, no. 2, pp. 1137–1167, 2019, doi: 

10.1007/s11071-019-05253-6. 

[15] S. Sharma and G. P. Samanta, “A Leslie-Gower predator-prey model with disease in 

prey incorporating a prey refuge,” Chaos, Solitons and Fractals, vol. 70, no. 1, pp. 

69–84, 2015, doi: 10.1016/j.chaos.2014.11.010. 

[16] S. Wang, Z. Ma, and W. Wang, “Dynamical behavior of a generalized eco-

epidemiological system with prey refuge,” Adv. Differ. Equations, vol. 2018, no. 1, 

pp. 1–20, 2018, doi: 10.1186/s13662-018-1704-x. 

[17] M. Onana, B. Mewoli, and J. J. Tewa, “Hopf bifurcation analysis in a delayed Leslie–

Gower predator–prey model incorporating additional food for predators, refuge and 

threshold harvesting of preys,” Nonlinear Dyn., vol. 100, no. 3, pp. 3007–3028, 2020, 

doi: 10.1007/s11071-020-05659-7. 

[18] Y. Cai, C. Zhao, W. Wang, and J. Wang, “Dynamics of a Leslie-Gower predator-

prey model with additive Allee effect,” Appl. Math. Model., vol. 39, no. 7, pp. 2092–

2106, 2015, doi: 10.1016/j.apm.2014.09.038.