## Analisis Kestabilan Model Predator-Prey dengan
Infeksi Penyakit pada Prey dan Pemanenan Proporsional pada Predator
Siti Maisaroh 1 , Resmawan 2 , Emli Rahmi 3,*
1,2,3 Jurusan Matematika, Universitas Negeri Gorontalo, Bone Bolango 96119, Indonesia
* Penulis Korespondensi: emlirahmi@ung.ac.id

## Jambura Journal of Biomathematics
Jambura J. Biomath. Volume 1, Issue 1, pp. 8-15, June 2020
Journal Homepage: http://ejurnal.ung.ac.id/index.php/JJBM/index https://doi.org/10.34312/jjbm.v1i1.5948
Received: 5 May 2020, Accepted: 16 June 2020, Published Online: 29 June 2020

## Abstrak
Artikel ini membahas dinamika dari model predator-prey dengan infeksi penyakit pada prey dan pemanenan pada predator . Prey
diasumsikan terbagi dalam dua kompartemen yaitu prey rentan dan prey terinfeksi. Model ini memiliki lima titik kesetimbangan yaitu titik kepunahan semua populasi, titik kepunahan prey terinfeksi dan predator , titik kepunahan prey terinfeksi, dan titik eksistensi semua populasi. Titik kepunahan semua populasi merupakan titik pelana sehingga kondisi ini tidak akan mungkin tercapai sedangkan untuk empat titik kesetimbangan lainnya stabil bersyarat. Diakhir artikel ini, untuk menunjang hasil analisis, diberikan simulasi numerik menggunakan metode Runge-Kutta Order empat.
Kata Kunci: Dinamika; Predator-prey ; Infeksi penyakit; Pemanenan

## Abstract
The dynamics of predator-prey model with infectious disease in prey and harvesting in predator is studied. Prey is divided into two compartments i.e the susceptible prey and the infected prey. This model has five equilibrium points namely the all population extinction point, the infected prey and predator extinction point, the infected prey extinction point, and the co-existence point. We show that all population extinction point is a saddle point and therefore this condition will never be attained, while the other equilibrium points are conditionally stable. In the end, to support analytical results, the numerical simulations are given by using the fourth-order Runge-Kutta method.
Keywords: Dynamics; Predator-prey; Infectious diseases; Harvesting

## 1. Pendahuluan
Pemodelan matematika merupakan cabang ilmu matematika yang merepresentasikan atau menjelaskan masalah di kehidupan nyata kedalam bentuk matematika. Pemodelan matematika biasanya selalu dikaitkan dengan cabang ilmu yang lain seperti biologi, fisika, kesehatan, dan teknik [1] . Salah satu cabang ilmu biologi adalah ekologi. Ekologi adalah ilmu yang mempelajari makhluk hidup dan interaksinya terhadap lingkungan maupun interaksi dengan sesama makhluk hidup [2] . Pada dasarnya model matematika dalam bidang ekologi merupakan studi tentang keterkaitan antara spesies dan lingkungannya, dalam bidang-bidang seperti interaksi pemangsa dan persaingan. Salah satu model matematika dalam bidang ekologi adalah model predator-prey . Model Lotka-Voltera adalah model predator-prey yang paling sederhana [3] . Berdasarkan model tersebut, dapat diketahui bahwa kedua spesies saling mempengaruhi satu sama lain, apabila terdapat spesies mangsa yang berlimpah, maka populasi pemangsa juga terus meningkat. Model Lotka-Voltera mengasumsikan laju pertumbuhan prey dan predator bertumbuh secara eksponensial. Kemudian pada tahun 1934, Gause mengembangkan model Lotka-Voltera [4] , pada model Gause tersebut laju pertumbuhan prey diasumsikan bertumbuh secara logistik dan pertumbuhan predator bertumbuh secara eksponensial [5] . Kajian tentang model predator-prey merupakan hal yang sangat menarik sehingga sampai saat ini masih terus dipelajari dan dikembangkan seperti di [6–13] .
Selain model ekologi, masalah epidemiologi merupakan salah satu topik hangat dalam pemodelan matematika. Pemodelan matematika dalam epidemiologi memberikan pemahaman tentang mekanisme mendasar yang mempengaruhi penyebaran penyakit. Dalam pemodelan matematika dari penularan penyakit, seperti di sebagian besar bidang pemodelan matematika lain, biasanya dilakukan modifikasi antara model-model sederhana, yang menghilangkan sebagian detail dan merancang untuk memodelkan perilaku kualitatif umum, dan model kompleks, biasanya dirancang untuk situasi tertentu [14] .
Pada beberapa tahun terakhir penelitian akan ekologi dan epidemiologi banyak dilakukan, meski dua bidang ini berbeda kontribusinya terhadap ilmu pengetahuan namun ada kesamaan diantara keduanya yaitu mencakup makhluk hidup. Model yang diangkat pada artikel ini merupakan gabungan kajian dari model Ekologi dan model Epidemiologi yang disebut dengan model Eko-epidemiologi. Beberapa peneliti telah mempelajari model eko-epidemiologi seperti pada [1, 15–21] . Artikel ini mempelajari model yang serupa dengan model dalam penelitian Panigoro dkk. [20] yaitu model eko-epidemiologi dengan mengasumsikan bahwa prey tumbuh logistik dan predator tumbuh secara eksponensial dengan adanya infeksi penyakit pada prey . Namun fungsi respon yang digunakan pada model penelitian [20] mengikuti Holling Tipe II sedangkan pada penelitian ini menggunakan fungsi respon Holling Tipe I. Selanjutnya pada beberapa kondisi yang terjadi di alam memperlihatkan terjadinya perburuan terhadap spesies predator oleh manusia. Oleh karena itu, model ini dimodifikasi dengan menambahkan asumsi bahwa terjadi pemanenan terhadap predator , yang berdasarkan kajian literatur penulis, model dengan modifikasi ini belum pernah dipelajari sebelumnya. Meskipun modifikasi yang dilakukan cukup sederhana, namun hasil analisis memperlihatkan bahwa terjadi perubahan signifikan pada dinamika dari model dibandingkan model rujukan. Eksistensi dan kestabilan dari model sangat bergantung pada besarnya pemanenan yang dilakukan. Hal inilah yang akan ditunjukkan dalam artikel ini.

## 2. Metode
Adapun langkah-langkah yang dilakukan dalam mempelajari model ini dijabarkan sebagai berikut.
1. Melakukan formulasi model berdasarkan asumsi dan kajian teoritis.
2. Mengidentifikasi titik kesetimbangan dan syarat kestabilannya.
3. Memberikan simulasi numerik yang bersesuaian dengan hasil analisis. Adapun skema numerik yang digunakan yaitu metode Runge-Kutta orde empat, yang dijalankan dengan aplikasi Python 3.7.

## 3. Hasil dan Pembahasan

## 3.1. Formulasi Model
Berikut adalah asumsi-asumsi yang digunakan dalam merumuskan model :
1. Dengan adanya penyakit pada prey menyebabkan prey terbagi atas dua kompartemen yaitu kompartemen rentan atau x S dan kompartemen terinfeksi atau x I , dimana prey terinfeksi tidak dapat kembali menjadi rentan.
2. Populasi prey rentan tumbuh logistik dengan laju pertumbuhan intrinsik r dan daya dukung lingkungan K . 3. Predator diasumsikan berinteraksi dengan prey rentan dan prey terinfeksi mengikuti fungsi respon Holling Tipe II.
4. Laju penyebaran penyakit antara prey rentan dan prey terinfeksi terjadi secara bilinear sebesar β . 5. Adanya kematian alami pada populasi prey terinfeksi dan predator .
6. Pemanenan terhadap predator terjadi secara proporsional sebesar h . 7. Tidak ada migrasi dalam populasi prey maupun predator .
Berdasarkan asumsi tersebut, dikonstruksi diagram kompartemen seperti pada Gambar 1 dan dirumuskan model seperti pada pers. (1) . Adapun interpretasi biologis dari setiap parameter dapat dilihat pada Tabel 1 .
dx S dt = rx S 1 − x S + x I K − m 1 x S y − β x S x I dx I dt = β x S x I − α x I − m 2 x I y dy dt = b 1 x S y + b 2 x I y − δ y − hy , (1)
3.2. Titik Kesetimbangan dan Kestabilannya
Didefinisikan R 3 + = { ( x S , x I , y ) | x S ≥ 0, x I ≥ 0, y ≥ 0, x S , x I , y ∈ R } . Titik kesetimbangan biologis E didapatkan dengan menyelesaikan dx S dt = dx I dt = dy dt = 0 dan memenuhi E ∈ R 3 + . Dengan demikian diperoleh lima titik kesetimbangan yaitu titik kepunahan semua populasi E 0 = ( 0, 0, 0 ) , titik kepunahan prey terinfeksi dan predator E 1 = ( K , 0, 0 ) , titik kepunahan prey terinfeksi E 2 = δ + h b 1 , 0, r 1 − δ + h b 1 K m 1 ! , titik kepunahan predator
Prey Rentan ( x S )
Prey Terinfeksi ( x I ) Predator ( y )
m 1 x s y m 2 x I y β x S x I rx S 1 − x S + x I K α x I δ y hy ( b 1 x S + b 2 x I ) y
Gambar 1. Diagram kompartemen model predator-prey dengan prey terinfeksi dan pemanenan pada predator
Tabel 1. Variabel dan parameter yang digunakan
Simbol Interpretasi Biologis x S ( t )
Prey rentan
x I ( t ) Prey terinfeksi
y ( t ) Predator r Laju pertumbuhan intrinsik prey rentan K Daya dukung lingkungan β Laju perpindahan kompartemen prey rentan ke prey terinfeksi m 1 Koefisien konversi prey rentan akibat adanya interaksi predator dengan prey rentan m 2 Koefisien konversi prey terinfeksi akibat adanya interaksi prey terinfeksi dengan predator
α Kematian alami pada prey terinfeksi b 1
Koefisien konversi pertumbuhan predator akibat adanya interaksi dengan prey rentan b 2
Koefisien konversi pertumbuhan predator akibat adanya interaksi dengan prey terinfeksi
δ Kematian alami predator h Pemanenan predator secara proporsional E 3 = α β , r ( β K − α ) β ( r + β K ) , 0 , titik eksistensi semua populasi E 4 = x ∗ S , δ + h − b 1 x
∗ S b 2 , β x ∗ S − α m 2 dengan
x ∗ S = rm 2 ( δ + h ) + Km 2 β ( δ + h ) − ( rKb 2 m 2 + Kb 2 m 2 α ) rb 1 m 2 + Km 2 b 1 − ( rb 2 m 2 + Kb 2 m 1 β ) .
Teorema 1. (i) E 2 eksis jika h < b 1 K − δ ; (ii) E 3 eksis jika α < β K; (iii) E 4 eksis jika x ∗ S > 0 dan α β < x ∗ S < δ + h b 1 .
bukti . Jika h < b 1 K − δ maka 1 − δ + h b 1 K > 0 sehingga E 2 ∈ R 3 + . Untuk E 3 , jika α < β K , maka β K − α > 0 sehingga E 3 ∈ R 3 + . Selanjutnya, jika x ∗ S > 0 dan α β < x ∗ S < δ + h b 1 mengakibatkan δ + h − b 1 x ∗ S > 0 dan β x ∗ S − α > 0 sehingga E 4 ∈ R 3 + .
Selanjutnya akan dipelajari dinamika disekitar titik-titik kesetimbangan tersebut. Untuk mengidentifikasi kestabilan lokal dari titik kesetimbangan, dilakukan pelinearan terhadap model (1) . Hasil peliniearan memberikan matriks Jacobi dari model (1) yaitu sebagai berikut.
J ( x S , x I , y ) =   r − 2 rx S K − rx I K − m 1 y − β x I rx S K − β x S − m 1 x S β x I β x S − α − m 2 y − m 2 x S b 1 y b 2 y b 1 x S + b 2 x I − δ − h   (2)
Syarat kestabilan dari suatu titik kesetimbangan diperoleh dengan mencari nilai eigen dari matriks Jacobi (2) yang dievaluasi pada setiap titik kesetimbangan.
Teorema 2. Titik kesetimbangan E 0 selalu merupakan titik pelana.
bukti . Dengan mensubtitusi E 0 ke matriks Jacobi (2) diperoleh
J ( E 0 ) =   r 0 0 0 − α 0 0 0 − ( δ + h )   , (3)
yang memberikan nilai eigen λ 1 = r , λ 2 = − α , dan λ 3 = − ( δ + h ) . Dengan demikian, E 0 merupakan titik kesetimbangan tidak stabil tipe pelana.
Teorema 3. Titik kesetimbangan E 1 stabil asimtotik lokal jika K < min n α β , δ + h b 1 o .
bukti . Titik kesetimbangan E 1 memiliki matriks Jacobi
J ( E 1 ) =   − r − K β − r − m 1 K 0 K β − α 0 0 0 Kb 1 − δ − h   , (4)
sehingga diperoleh nilai eigen λ 1 = − r , λ 2 = K β − α , dan λ 3 = Kb 1 − δ − h . Jika K < min n α β , δ + h b 1 o maka K β − α < 0 dan Kb 1 − δ − h yang mengakibatkan λ 1,2,3 < 0. Akibatnya E 1 stabil asimtotik lokal.
Teorema 4. Misalkan
η 1 = K α b 1 + Krb 1 m 2 , η 2 = K βδ m 1 + K β m 1 h + δ rm 2 + hrm 2 , D 1 = r ( δ + h ) Kb 1 2 − 4 r ( δ + h )( Kb 1 − ( δ + h )) Kb 1 .
Titik kesetimbangan E 2 stabil asimtotik lokal jika η 1 > η 2 dan
(i) D 1 < 0 , atau; (ii) D 1 ≥ 0 dan K > δ + h b 1 .
bukti . Matriks Jacobi model (1) pada titik E 2 adalah
J ( E 2 ) =     − r ( δ + h ) b 1 k − ( δ + h )( K β + r ) Kb 1 − m 1 ( δ + h ) b 1 0 − η 1 − η 2 Km 1 b 1 0 r ( Kb 1 − δ − h ) Km 1 r ( Kb 1 − δ − h ) b 2 Km 1 b 1 0     (5)
sehingga diperoleh nilai eigen λ 1 = − η 1 − η 2 Km 1 b 1 , λ 2 = 1 2 − r ( δ + h ) b 1 K + √ D 1 , dan λ 3 = 1 2 − r ( δ + h ) b 1 K − √ D 1 . Jika η 1 > η 2 maka λ 1 < 0 sehingga kestabilannya bergantung pada nilai λ 2,3 . Karena Re ( λ 2,3 ) < 0 maka kestabilannya bergantung pada D 1 . Jika D 1 < 0 maka λ 2,3 ∈ C dengan Re ( λ 2,3 ) < 0. Dengan demikian, E 2 stabil asimtotik lokal. Jika D 1 ≥ 0 maka λ 2,3 ∈ R . Agar λ 2,3 < 0, dengan mudah dapat ditunjukkan bahwa K > δ + h b 1 harus dipenuhi. Dengan demikian E 2 stabil asimtotik lokal.
Teorema 5. Misalkan
ξ 1 = K α b 1 β + α rb 1 + K β b 2 r ,
ξ 2 = K β 2 δ + K β 2 rb 2 + α rb 2 + βδ r + β hr , D 2 = r α β K 2 − 4 α r ( K β + r ( K β − α )) K β ( K β + r ) .
Titik kesetimbangan E 3 stabil asimtotik lokal jika ξ 1 > ξ 2 dan; (i) D 2 < 0 , atau; (ii) D 2 ≥ 0 dan K > α β .
x S 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 x I 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 y 0.0 0.1 0.2 0.3 0.4 0.5 E 0
Gambar 2. Simulasi numerik disekitar titik kesetimbangan E 0 dengan nilai parameter (8)
bukti . Titik kesetimbangan E 3 memiliki matriks Jacobi
J ( E 3 ) =     − r α β K − α ( K β + r ) K β − m 1 α β r ( K β − α ) K β + r 0 − r ( K β − α ) m 2 β ( K β + r )
0 0 ξ 1 − ξ 2 β ( K β + r )     , (6)
sehingga didapatkan nilai eigen λ 1 = ξ 1 − ξ 2 β ( K β + r ) , λ 2 = − 1 2 r α β K − √ D 2 , dan λ 3 = − 1 2 r α β K + √ D 2 . Jika ξ 1 < ξ 2 maka λ 1 < 0. Perhatikan bahwa Re ( λ 2,3 ) < 0 sehingga kestabilan E 3 bergantung pada nilai D 2 . Jika D 2 < 0 maka λ 2,3 ∈ C dengan Re ( λ 2 , 3 ) < 0. Oleh karena itu E 3 stabil asimtotik lokal. Jika D 2 ≥ 0 maka λ 2,3 ∈ R . Jika K > α β maka λ 2,3 < 0. Dengan demikian syarat kestabilan untuk E 3 telah dibuktikan.
Teorema 6. Misalkan
K 1 = m 1 b 2 x ∗ S ( β x ∗ S − α ) m 2 + β x ∗ S ( r + β K )( δ + h − b 1 x ∗ S ) Kb 2 + ( δ + h − b 1 x ∗ S )( β x
∗ S − α )
K 2 = m 1 β x ∗ S ( β x ∗ S − α )( δ + h − b 1 x ∗ S ) m 2 + rx ∗ S ( β x ∗ S − α )( δ + h − b 1 x ∗ S ) K − b 1 x ∗ S ( r + β K )( δ + h − b 1 x ∗ S )( β x ∗ S − α ) Kb 2
Titik kesetimbangan E 3 stabil asimtotik lokal jika 0 < K 2 < rx ∗ S K 1 K .
bukti . Matriks jacobi untuk titik kesetimbangan E 4 = x ∗ S , δ + h − b 1 x
∗ S b 2 , β x ∗ S − α m 2 yaitu
J ( E 4 ) =     − rx ∗ S K − x ∗ S K ( r + β K ) − m 1 x ∗ S β ( δ + h − b 1 x ∗ S ) b 2 0 − m 2 ( δ + h − b 1 x ∗ S ) b 2 b 1 ( β x ∗ S − α ) m 2 b 2 ( β x ∗ S − α ) m 2 0     , (7)
yang memberikan polinom karakteristik P ( λ ) = λ 3 + rx ∗ S K λ 2
+ K 1 λ + K 2 . Dengan menggunakan kriteria kestabilan Routh-Hurwitz, diperoleh syarat kestabilan E 4 yaitu 0 < K 2 < rx ∗ S K 1 K .

## 3.3. Simulasi Numerik
Untuk memperkuat hasil analisis yang diberikan pada sebelumnya, dilakukan simulasi numerik terhadap model (1) dengan menggunakan metode Runge-Kutta orde empat. Karena keterbatasan data yang ada, pemilihan parameter disesuaikan dengan kondisi kestabilan yang diberikan pada hasil sebelumnya. Dipilih parameter awal sebagai berikut.
r = 0.4, m 1 = 0.3, β = 0.1, m 2 = 0.2, α = 0.2, b 1 = 0.2, b 2 = 0.2, δ = 0.1, h = 0.2, K = 1. (8)
x S 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 x I 0.00 0.02 0.04 0.06 0.08 0.10 y 0.0 0.1 0.2 0.3 0.4 0.5 E 1 (a) 0 10 20 30 40 50 t 0.0 0.2 0.4 0.6 0.8 1.0 Populasi x S x I y (b)
Gambar 3. Simulasi numerik disekitar titik kesetimbangan E 1 dengan nilai parameter (8)
x S 0.40 0.45 0.50 0.55 0.60 0.65 0.70 x I 0.000 0.025 0.050 0.075 0.100 0.125 0.150 0.175 0.200 y 0.40 0.45 0.50 0.55 0.60 0.65 0.70 E 2 (a) 0 10 20 30 40 50 60 70 t 0.0 0.1 0.2 0.3 0.4 0.5 0.6 Populasi x S x I y (b)
Gambar 4. Simulasi numerik disekitar titik kesetimbangan E 2 dengan nilai parameter (8) dan h = 0.01
Berdasarkan Teorema 2 , titik kesetimbangan E 0 selalu tidak stabil tipe pelana seperti pada Gambar 2 . Dengan nilai parameter yang sama seperti pada persamaan (8) , berdasrkan Teorema 3 , terdapat titik kesetimbangan E 1 yang stabil asimtotik lokal, lihat Gambar 3 . Populasi prey terinfeksi dan predator menuju kepunahan, sehingga yang tersisa hanyalah populasi prey yang rentan.
Selanjutnya diberikan simulasi numerik disekitar titik kesetimbangan E 2 . Agar kondisi kestabilan seperti pada Teorema 4 terpenuhi, nilai dari h pada parameter (8) diturunkan menjadi 0.01. Akibatnya semua solusi yang cukup dekat dengan titik kesetimbangan E 2 akan menuju ke titik kesetimbangan tersebut, seperti pada Gambar 4 . Hal ini berarti populasi prey terinfeksi akan punah, sedangkan eksistensi prey rentan dan predator akan terjaga.
Pada saat β naik menjadi 0.3, muncul titik kesetimbangan E 3 yang stabil asimtotik lokal. Berdasarkan Teorema 5 , untuk nilai awal yang cukup dekat dengan E 3 akan konvergen ke E 3 seperti yang ditunjukkan pada Gambar 5 . Populasi predator akan punah sedangkan populasi prey terjaga, namun infeksi penyakit tetap terjadi.
Di akhir simulasi, akan ditunjukkan kestabilan dari titik eksistensi semua populasi E 4 . Dipilih nilai parameter sebagai berikut.
r = 0.3, m 1 = 3.03, β = 0.05, m 2 = 1, α = 0.2, b 1 = 0.5, b 2 = 2, δ = 2, h = 0.63, K = 14.5. (9)
Berdasarkan Teorema 6 , solusi berosilasi menuju titik kesetimbangan E 4 . Hal ini berarti solusi dengan nilai awal yang cukup dekat dengan E 4 mengakibatkan eksistensi dari seluruh populasi terjaga, lihat Gambar 6 . Dari kelima gambar diatas, satu-satunya kompartemen yang tidak mungkin punah adalah prey rentan, sedangkan lainnya dapat punah dengan kondisi parameter tertentu.
x S 0.60 0.62 0.64 0.66 0.68 0.70 0.72 x I 0.170 0.175 0.180 0.185 0.190 0.195 0.200 0.205 0.210 y 0.000 0.002 0.004 0.006 0.008 0.010 0.012 E 3 (a) 0 50 100 150 200 250 300 t 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 Populasi x S x I y (b)
Gambar 5. Simulasi numerik disekitar titik kesetimbangan E 3 dengan nilai parameter (8) , h = 0.01, dan β = 0.3
x S 5.2612 5.2614 5.2616 5.2618 5.2620 x I 1e 8+3.718e 4 4.0 4.5 5.0 5.5 6.0 6.5 7.0 y 0.063060 0.063065 0.063070 0.063075 0.063080 0.063085 0.063090 0.063095 E 3 (a) 0 50 100 150 200 250 300 t 5.26125 5.26150 5.26175 5.26200 x S 0 50 100 150 200 250 300 t 4 5 6 7 x I 1e 8+3.718e 4 0 50 100 150 200 250 300 t 0.06306 0.06307 0.06308 0.06309 y (b)
Gambar 6. Simulasi numerik disekitar titik kesetimbangan E 3 dengan nilai parameter (9)

## 4. Kesimpulan
Dinamika model predator-prey dengan infeksi penyakit pada prey dan pemanenan pada predator telah dipelajari. Hasil analisis menunjukkan bahwa terdapat empat titik kesetimbangan yaitu titik kepunahan semua populasi, titik kepunahan prey terinfeksi dan predator , titik kepunahan prey terinfeksi, dan titik eksistensi semua populasi. Telah ditunjukkan bahwa kepunahan seluruh populasi tidak akan mungkin tercapai, selain itu juga diperlihatkan bahwa prey rentan adalah satu-satunya populasi yang tidak akan mungkin punah pada model ini. Secara analitik telah ditunjukkan kondisi yang harus dipenuhi agar setiap titik kesetimbangan stabil asimtotik lokal. Dari hasil analisis juga telah ditunjukkan adanya kondisi yang dapat dipenuhi agar populasi predator dan prey dapat dipertahankan, dengan infeksi penyakit pada prey dapat dihilangkan. Diakhir artikel telah ditunjukkan simulasi numerik berupa potret fase dan deret waktu disekitar setiap titik kesetimbangan serta interpretasi biologisnya.

## Ucapan Terimakasih
Terimakasih kami ucapkan kepada Bapak Hasan S. Panigoro yang telah banyak membantu dalam penyelesaian artikel ini.

## Referensi
[1] S. Kant and V. Kumar, “Stability analysis of predator–prey system with migrating prey and disease infection in both species,” Applied Mathematical Modelling , vol. 42, pp. 509–539, 2017.
[2] J. B. Reece and N. A. Campbell, Campbell Biology , 9th ed. Pearson Australia, 2011.
[3] A. J. Lotka, “Elements of Physical Biology,” Nature , vol. 116, no. 2917, pp. 461–461, 1925.
[4] A. A. Berryman, “The Orgins and Evolution of Predator-Prey Theory,” Ecology , vol. 73, no. 5, pp. 1530–1535, 1992.
[5] J. D. Murray, Mathematical Biology I: An Introduction , 3rd ed., ser. Interdisciplinary Applied Mathematics, S. Antman, J. Marsden, L. Sirovich, and S. Wiggins, Eds. New York, NY: Springer New York, 2002.
[6] H. S. Panigoro and E. Rahmi, “Modifikasi sistem predator-prey: dinamika model Leslie-Gower dengan daya dukung yang tumbuh logistik,” in SEMIRATA MIPAnet . UNSRAT Manado, 2017, pp. 94–103.
[7] F. M. Hilker and E. Liz, “Proportional threshold harvesting in discrete-time population models,” Journal of Mathematical Biology , vol. 79, no. 5, pp. 1927–1951, 2019.
[8] A. Suryanto, I. Darti, H. S. Panigoro, and A. Kilicman, “A Fractional-Order Predator–Prey Model with Ratio-Dependent Functional Response and Linear Harvesting,” Mathematics , vol. 7, no. 11, p. 1100, 2019.
[9] N. Hasan, R. Resmawan, and E. Rahmi, “Analisis Kestabilan Model Eko-Epidemiologi dengan Pemanenan Konstan pada Predator,” Jurnal Matematika, Statistika dan Komputasi , vol. 16, no. 2, p. 121, 2019.
[10] D. Mukherjee, “Study of fear mechanism in predator-prey system in the presence of competitor for the prey,” Ecological Genetics and Genomics , vol. 15, p. 100052, 2020.
[11] J. Alidousti and E. Ghafari, “Dynamic behavior of a fractional order prey-predator model with group defense,” Chaos, Solitons and Fractals , vol. 134, p. 109688, 2020.
[12] H. S. Panigoro, E. Rahmi, N. Achmad, and S. L. Mahmud, “The Influence of Additive Allee Effect and Periodic Harvesting to the Dynamics of Leslie-Gower Predator-Prey Model,” Jambura Journal of Mathematics , vol. 2, no. 2, pp. 87–96, 2020.
[13] S. H. Arsyad, R. Resmawan, and N. Achmad, “Analisis Model Predator-Prey Leslie-Gower dengan Pemberian Racun Pada Predator,” Jurnal Riset dan Aplikasi Matematika (JRAM) , vol. 4, no. 1, pp. 1–16, 2020.
[14] F. Brauer and C. Castillo-Chavez, Mathematical Models in Population Biology and Epidemiology , ser. Texts in Applied Mathematics. Springer New York, 2001.
[15] S. Rida, M. Khailil, H. A. Hosham, and S. Gadellah, “Predator-Prey Fractional-Order Dynamical System with Both the Populations Affected by Diseases,” Journal of Fractal Calculus and Applications , vol. 5, no. 13, pp. 1–11, 2014.
[16] A. Suryanto, I. Darti, and S. Anam, “Stability analysis of pest-predator interaction model with infectious disease in prey,” 2018, p. 020018.
[17] A. Mondal, A. K. Pal, and G. P. Samanta, “On the dynamics of evolutionary Leslie-Gower predator-prey eco-epidemiological model with disease in predator,” Ecological Genetics and Genomics , vol. 10, no. August 2018, p. 100034, 2019.
[18] A. S. Purnomo, I. Darti, and A. Suryanto, “Dynamics of eco-epidemiological model with harvesting,” AIP Conference Proceedings , vol. 1913, no. December 2017, 2017.
[19] A. Sha, S. Samanta, M. Martcheva, and J. Chattopadhyay, “Backward bifurcation, oscillations and chaos in an eco-epidemiological model with fear effect,” Journal of biological dynamics , vol. 13, no. 1, pp. 301–327, 2019.
[20] H. S. Panigoro, A. Suryanto, W. M. Kusumawinahyu, and I. Darti, “Dynamics of a Fractional-Order Predator-Prey Model with Infectious Diseases in Prey,” Commun. Biomath. Sci. , vol. 2, no. 2, pp. 105–117, 2019.
[21] M. Moustafa, M. H. Mohd, A. I. Ismail, and F. A. Abdullah, “Dynamical analysis of a fractional-order eco- epidemiological model with disease in prey population,” Advances in Difference Equations , vol. 2020, no. 1, p. 48,
2020.
©  2020  by  the Authors .  This  article  is  an  open  access  article  distributed  under  the  terms  and  conditions  of  the  Creative Commons  A � ribution-NonComercial  4.0  International  License. Editorial  of  JJBM: Department  of  Mathematics,  State University  of  Gorontalo,  Jln.  Prof.  Dr.  Ing.  B.  J.  Habibie,  Moutong,  Tilongkabila,  Kabupaten  Bone  Bolango,  Provinsi Gorontalo 96119, Indonesia.