## PENYELESAIAN MODIFIKASI MODEL PREDATOR PREY LESLIE-GOWER DENGAN SEBAGIAN PREY

## TERINFEKSI MENGGUNAKAN ADAMS BASHFORTH MOULTON ORDE EMPAT
( Solving of Modified Leslie-Gower Predator-Prey Model with Prey Infection Using Fourth Order Adams Bashforth Moulton )

## Liatri Arianti, Rusli Hidayat, Kosala Dwidja Purnomo
Jurusan Matematika, Fakultas MIPA, Universitas Jember Jl. Kalimantan 37 Jember 68121, Indonesia E-mail: {rusli, kosala}.fmipa@unej.ac.id
Abstract. Eco-epidemiology  is  a  science  that  studies  the  spread  of  infectious diseases in a population in an ecosystem where two or more species interact like a predator prey. In this paper discusses about how to solve modification Leslie Gower of predator prey models (with Holling II response function) with some prey infected using  fourth  order  Adams  Bashforth  Moulton  method.  This  paper  used  a  simple disease-spreading model that is Susceptible-Infected (SI). The model is divided into three  populations:  the  sound  prey  (which  is  susceptible),  the  infected  prey  and predator population.
Keywords: Adams Basforth Moulton, Eco-epidemiology Holling Tipe II, Local stability, Leslie-Gower, Predator-Prey model
MSC 2010: 65C20

## 1. Pendahuluan
Salah satu model interaksi antar mahkluk hidup adalah predator prey. Alfred dan Vito Volterra pada pertengahan tahun 1920 mengembangkan sepasang persamaan diferensial sederhana  yang  menggambarkan  fenomena predator  prey untuk  pertama  kali  dikenal sebagai model Lotka-Volterra [1]. Model sederhana tersebut kemudian dimodifikasi oleh Leslie  dan  Gower  dengan  memperkenalkan  model predator  prey dimana carrying capacity lingkungan predator sebanding dengan banyaknya prey [2]. Holling pada tahun 1959 memperkenalkan fungsi respons. Holling membagi fungsi respons menjadi tiga tipe, yaitu fungsi respons tipe I, tipe II, dan tipe III [3].   Masalah lain yang dapat terjadi pada suatu populasi adalah terjadinya epidemi. Chattopadhyay dkk [4] telah mengembangkan model  matematika predator  prey dengan  mangsa  yang  terinfeksi.  Zhou  dkk  [5]  telah mengusulkan  modifikasi  model predator  prey Leslie  Gower  (dengan  fungsi  respon Holling II) dengan sebagian prey terinfeksi. Zhou dkk menggunakan model penyebaran
penyakit  yang  sederhana  yaitu Susceptible-Infected (SI).  Zhou  dkk  membagi prey ke dalam dua kelompok, yaitu kelompok rentan dan kelompok terinfeksi.
Metode  Adams  Bashforth  Moulton  adalah  metode  penyelesaian  diferensial  banyak langkah karena untuk menyelesaikannya diperlukan lebih dari satu informasi nilai awal, sehingga metode ini harus menggunakan metode lain untuk mencari nilai awal. Metode Runge-Kutta  dipilih  untuk  pengambilan  nilai  awal  karena  pada  penelitian  sebelumnya oleh Putri [6] dalam skripsinya telah menyelesaikan model predator prey klasik dengan menggunakan  metode  Runge-Kutta  orde  empat  dan  Gill.  Penelitian  tersebut menyimpulkan  kedua  metode  tersebut  baik  dalam  penyelesain  model.  Metode  Adams Bashforth  Moulton  Orde  Empat  adalah  metode  yang  banyak  digunakan  karena mempunyai  ketelitian  yang tinggi.  Mardiana [7] dalam skripsinya  yang  menyelesaikan model predator prey Holling dengan faktor pemanenan, Mardiana menyimpulkan bahwa hasil  yang  diperoleh  dari  penyelesaian  numerik  dengan  menggunakan  metode  Adams Bashforth  Moulton  orde  empat  sama  dengan  hasil  analitik  pada  titik  kesetimbangan. Penulis tertarik untuk menganalis kestabilan model  yang di usulkan oleh Zhou dkk [5] dan  menyelesaikannya  dengan  menggunakan  metode  Adams  Bashforth  Moulton  orde empat.
Modifikasi Model Predator Prey Leslie Gower dengan Sebagian Prey Terinfeksi Pemodelan matematika epidemi telah menjadi topik penelitian yang menarik dan penting. Pada model ini akan dipelajari pengaruh penyakit pada lingkungan di mana terdapat dua atau  lebih  spesies  yang  saling  berinteraksi.  Model  dalam  tulisan  ini  terdiri  dari  tiga populasi, yaitu populasi prey rentan, prey terinfeksi, dan predator. Penyakit pada model ini hanya menyebar di antara populasi prey saja dan penyakit ini tidak diwariskan secara genetis. Populasi yang terinfeksi tidak pulih atau menjadi kebal. Populasi predator hanya memakan prey yang terinfeksi dengan skema Leslie-Gower Holling-type II.
                                2 2 2 1 1 1 K I Iy c a y K I Iy c cI SI I SI K I S rS S     
(1)
dengan kondisi awal 𝑆(0) = 𝑆 0 𝑆(𝑡) = kepadatan populasi prey rentan pada saat t 𝐼(𝑡) = kepadatan populasi prey terinfeksi pada saat t 𝑦(𝑡) = kepadatan populasi predator pada saat
𝑟 = laju pertumbuhan intrinsik populasi prey
𝐾 = kapasitas daya tampung prey rentan
𝛽 = laju infeksi (penularan penyakit) 𝑐 = laju kematian prey terinfeksi 𝑎 2 = laju pertumbuhan predator
𝑐 1 = nilai maksimum laju per kapita pengurangan I karena y 𝑐 2 = nilai maksimum laju per kapita pengurangan y karena I 𝐾 1 = daya dukung lingkungan untuk prey terinfeksi 𝐾 2 = daya dukung lingkungan untuk predator dan semua parameter bernilai positif [5].

## 2. Metodologi
Metode penelitian yang digunakan untuk penyelesaian modifikasi model Leslie-Gower dengan  sebagian prey terinfeksi  menggunakan Adams  Bashorth  Moulton  orde  empat adalah sebagai berikut:
1. Menentukan titik kesetimbangan model dengan menyelesaikan Persamaan (1) sama dengan nol.
𝑓 = 𝑆̇ = 𝑟𝑆(1 − 𝑆+𝐼 𝐾 ) − 𝛽𝑆𝐼 = 0 𝑔 = 𝐼̇ = 𝛽𝑆𝐼 − 𝑐𝐼 − 𝑐 1 𝐼𝑦 𝐼+𝐾 1 = 0 ℎ = 𝑦̇ = (𝑎 2 − 𝑐 1 𝐼𝑦 𝐼+𝐾 2 )𝑦 = 0
dengan 𝑆 * , 𝐼 * , 𝑦 * merupakan nilai kesetimbangan dari model.
2. Linierisasi model menggunakan matriks Jacobian, kemudian menentukan kestabilan di sekitatar titik kesetimbangan. Matriks Jacobian dari Persamaan (1) adalah sebagai berikut:
𝐽 = ( 𝑎 11 𝑎 12 0 𝛽𝐼 𝑎 22 𝑎 23 0 𝑎 32 𝑎 33 ) (2)
dengan
𝑎 11 = 𝑟 (1 − 𝑆 + 𝐼 𝐾 ) − 𝑟𝑆 𝐾 − 𝛽𝐼 𝑎 12 = −𝛽𝑆 − 𝑟𝑆 𝐾 𝑎 22 = 𝛽𝑆 − 𝑐 − 𝑐 1 𝑦 𝐼 + 𝐾 1 + 𝑐 1 𝐾 1 𝑦 (𝐼 + 𝐾 1 ) 2 𝑎 23 = −𝑐 1 𝐼 (𝐼 + 𝐾 1 ) 𝑎 32 = 𝑐 2 𝑦 2 (𝐼 + 𝐾 2 ) 2 𝑎 33 = 𝑎 2 − 2𝑐 2 𝑦 𝐼 * + 𝐾 2
3. Penyelesaian  model  menggunakan  metode  numerik  Adams  Bashorth  Moulton  orde empat,  kemudian  membuat  program  menggunakan software MATLAB.  R2015b. Pengambilan  nilai  awal  menggunakan  Runge  Kutta  orde  4,  dengan  rumus  sebagai berikut:
𝑥 𝑛+1 = 𝑥 𝑛 + 1 6 (𝑘 1 + 2𝑘 2 + 2𝑘 3 + 𝑘 4 )ℎ (3)
dengan,
𝑘 1 = 𝑓(𝑡 𝑖 , 𝑥 𝑖 ) 𝑘 2 = 𝑓(𝑡 𝑖 + 1 2 ℎ, 𝑥 𝑖 + 1 2 ℎ𝑘 1 ) 𝑘 3 = 𝑓 (𝑡 𝑖 + 1 2 ℎ, 𝑥 𝑖 + 1 2 ℎ𝑘 2 ) 𝑘 4 = 𝑓(𝑡 𝑖 + ℎ, 𝑥 𝑖 + ℎ𝑘 3 )
Dilanjutkan dengan menggunakan Adams Bashforth Moulton orde empat.
Predictor :
𝑦 𝑖+1 = 𝑦 𝑖 + ℎ 24 (−9𝑓 𝑖−3 + 37𝑓 𝑖−2 − 59𝑓 𝑖−1 + 55𝑓 𝑖 ) (4)
Corrector:
𝑦 * 𝑖+1 = 𝑦 𝑖 + ℎ 24 (𝑓 𝑖−2 − 5𝑓 𝑖−1 + 19𝑓 𝑖 + 9𝑓 𝑖+1 ) (5)

## 3. Hasil dan Pembahasan
Dari penelitian ini didapat hasil sebagai berikut:

## Titik Kesetimbangan dan Kestabilannya
1. Titik Setimbang Kepunahan Setiap Populasi yaitu 𝑆 * = 𝐼 * = 𝑦 * = 0 didapatkan 𝐸 1 = (0,0,0) bersifat tidak stabil.
2. Titik Setimbang Prey Rentan dan Prey Terinfeksi Punah yaitu 𝑆 * = 𝐼 * = 0 didapatkan 𝐸 2 = (0,0, 𝑎 2 𝐾 2 𝑐 2 ) bersifat tidak stabil.
3. Titik Setimbang Prey Terinfeksi  dan Predator Punah  yaitu 𝐼 * = 𝑦 * = 0 didapatkan 𝐸 3 = (𝐾, 0,0) bersifat tidak stabil.
4. Titik Setimbang Prey Terinfeksi Mengalami Kepunahan yaitu 𝐼 * = 0 didapatkan 𝐸 4 = (𝐾, 0, 𝑎 2 𝐾 2 𝑐 2 ) bersifat stabil asimtotik dengan syarat 𝛽 < 𝑐𝑐 2 𝐾 1 +𝑐 1 𝑎 2 𝐾 2 𝐾𝑐 2 𝐾 1 .
5. Titik  Setimbang Predator Punah  yaitu 𝑦 * = 0 didapat 𝐸 5 = ( 𝑐 𝛽 , 𝑟(𝛽𝐾−𝑐) 𝛽(𝑟−𝛽𝐾) , 0) bersifat tidak stabil.
6. Titik Setimbang Tidak Ada Kepunahan didapat 𝐸 * = (𝑆 * , 𝐼 * , 𝑦 * ) yaitu
𝐼 + * (𝑐 2 𝐾 1 (𝛽 + 𝛽 2 𝐾 𝑟 ) + 𝑐 1 𝑎 2 − 𝑐 2 (𝛽𝐾 − 𝑐)) + √∆
2𝑐 2 (𝛽 + 𝛽 2 𝐾 𝑐 ) 𝑦 * = 𝑎 2 (𝐼 * + 𝐾 2 ) 𝑐 2 𝑆 * = 1 𝛽 (𝑐 + 𝑐 2 𝑦 * 𝐼 * + 𝐾 1 )
Bersifat stabil asimtotik dengan syarat
𝑐 2 𝐼 * +𝐾 2 > 𝑐 1 𝐼 * (𝐼 * +𝐾 1 ) 2 dan 𝐾 1 > 𝐾 2 .

## Penyelesaian Model dengan Adams Bashforth Moulton Orde Empat
Pencarian nilai awal menggunakan metode Runge Kutta orde empat pada Persamaan (1) menggunakan rumus Persamaan (3):
𝑆 𝑛+1 = 𝑆 𝑛 + 1 6 (𝑘 1 + 2𝑘 2 + 2𝑘 3 + 𝑘 4 )ℎ
𝐼 𝑛+1 = 𝐼 𝑛 + 1 6 (𝑙 1 + 2𝑙 2 + 2𝑙 3 + 𝑙 4 )ℎ
𝑦 𝑛+1 = 𝑦 𝑛 + 1 6 (𝑜 1 + 2𝑜 2 + 2𝑜 3 + 𝑜 4 )ℎ
Persamaan diferensial pertama
𝑘 1 = 𝑓(𝑡 1 , 𝑆 1 , 𝐼 1 , 𝑦 1 )
𝑘 2 = 𝑓 (𝑡 2 + ℎ 2 , 𝑆 2 + ℎ 2 𝑘 1 , 𝐼 2 + ℎ 2 𝑙 1 , 𝑦 2 + ℎ 2 𝑜 1 ) 𝑘 3 = 𝑓 (𝑡 3 + ℎ 2 , 𝑆 3 + ℎ 2 𝑘 2 , 𝐼 2 + ℎ 2 𝑙 2 , 𝑦 2 + ℎ 2 𝑜 2 ) 𝑘 4 = 𝑓(𝑡 4 + ℎ, 𝑆 4 + ℎ𝑘 3 , 𝐼 2 + ℎ𝑙 3 , 𝑦 2 + ℎ𝑜 3 )
Persamaan diferensial kedua
𝑙 1 = 𝑔(𝑡 1 , 𝑆 1 , 𝐼 1 , 𝑦 1 )
𝑙 2 = 𝑔 (𝑡 2 + ℎ 2 , 𝑆 2 + ℎ 2 𝑘 1 , 𝐼 2 + ℎ 2 𝑙 1 , 𝑦 2 + ℎ 2 𝑜 1 ) 𝑙 3 = 𝑔 (𝑡 3 + ℎ 2 , 𝑆 3 + ℎ 2 𝑘 2 , 𝐼 2 + ℎ 2 𝑙 2 , 𝑦 2 + ℎ 2 𝑜 2 ) 𝑙 4 = 𝑔(𝑡 4 + ℎ, 𝑆 4 + ℎ𝑘 3 , 𝐼 2 + ℎ𝑙 3 , 𝑦 2 + ℎ𝑜 3 )
Persamaan diferensial ketiga
𝑜 1 = ℎ(𝑡 1 , 𝑆 1 , 𝐼 1 , 𝑦 1 )
𝑜 2 = ℎ (𝑡 2 + ℎ 2 , 𝑆 2 + ℎ 2 𝑘 1 , 𝐼 2 + ℎ 2 𝑙 1 , 𝑦 2 + ℎ 2 𝑜 1 ) 𝑜 3 = ℎ (𝑡 3 + ℎ 2 , 𝑆 3 + ℎ 2 𝑘 2 , 𝐼 2 + ℎ 2 𝑙 2 , 𝑦 2 + ℎ 2 𝑜 2 ) 𝑜 4 = ℎ(𝑡 4 + ℎ, 𝑆 4 + ℎ𝑘 3 , 𝐼 2 + ℎ𝑙 3 , 𝑦 2 + ℎ𝑜 3 )
setelah didapatkan nilai-nilai awal yang dibutuhkan, selanjutnya dilakukan tahap prediksi menggunakan  skema  Adams  Basforth  Moulton  orde  empat  yaitu  Persamaan  (4)  pada Persamaan (1). Sehingga diperoleh:
Persamaan diferensial satu
𝑆 𝑟+1 = 𝑆 𝑟 + ℎ 24 (−9f(𝑡 𝑟−3 , 𝑆 𝑟−3 , 𝐼 𝑟−3 , 𝑦 𝑟−3 )) +37f(𝑡 𝑟−2 , 𝑆 𝑟−2 , 𝐼 𝑟−2 , 𝑦 𝑟−2 ) −59f(𝑡 𝑟−1 , 𝑆 𝑟−1 , 𝐼 𝑟−1 , 𝑦 𝑟−1 ) +55f(𝑡 𝑟 , 𝑆 𝑟 , 𝐼 𝑟 , 𝑦 𝑟 )
Persamaan diferensial dua
𝐼 𝑟+1 = 𝐼 𝑟 + ℎ 24 (−9g(𝑡 𝑟−3 , 𝑆 𝑟−3 , 𝐼 𝑟−3 , 𝑦 𝑟−3 )) +37g(𝑡 𝑟−2 , 𝑆 𝑟−2 , 𝐼 𝑟−2 , 𝑦 𝑟−2 ) −59g(𝑡 𝑟−1 , 𝑆 𝑟−1 , 𝐼 𝑟−1 , 𝑦 𝑟−1 ) +55g(𝑡 𝑟 , 𝑆 𝑟 , 𝐼 𝑟 , 𝑦 𝑟 )
persamaan differential tiga
𝑦 𝑟+1 = 𝑦 𝑟 + ℎ 24 (−9h(𝑡 𝑟−3 , 𝑆 𝑟−3 , 𝐼 𝑟−3 , 𝑦 𝑟−3 )) +37h(𝑡 𝑟−2 , 𝑆 𝑟−2 , 𝐼 𝑟−2 , 𝑦 𝑟−2 )
−59h(𝑡 𝑟−1 , 𝑆 𝑟−1 , 𝐼 𝑟−1 , 𝑦 𝑟−1 ) +55h(𝑡 𝑟 , 𝑆 𝑟 , 𝐼 𝑟 , 𝑦 𝑟 )
Pada  tahap  selanjutnya  akan  dilakukan  tahap  koreksi  dengan  menggunakan  Adams Bashforth Moulton orde empat pada Persamaan (5). Sehingga diperoleh:
Persamaan diferensial satu
𝑆 * 𝑟+1 = 𝑆 𝑟 + ℎ 24 (𝑓 𝑟−2 − 5𝑓 𝑟−1 + 19𝑓 𝑟 + 9𝑓 𝑟+1 )
Persamaan diferensial dua
𝐼 * 𝑟+1 = 𝐼 𝑟 + ℎ 24 (𝑔 𝑟−2 − 5𝑔 𝑟−1 + 19𝑔 𝑟 + 9𝑔 𝑟+1 )
Persamaan diferensial tiga
𝑦 * 𝑟+1 = 𝑦 𝑟 + ℎ 24 (ℎ 𝑟−2 − 5ℎ 𝑟−1 + 19ℎ 𝑟 + 9ℎ 𝑟+1 )
Hasil  diskritisasi  diatas  kemudian  di  implementasikan  dalam  pembuatan  program menggunakan software MATLAB R2015b.

## Simulasi dan Interpretasi
Nilai  parameter  yang  digunakan  dalam  simulasi  ini  sesuai  dalam  jurnal  Zhou  dkk  [5] dengan  nilai  awal 𝑆(0) = 1, 𝐼(0) = 0.5, 𝑦(0) = 1.2 dengan  waktu 𝑡 = 0 sampai 𝑡 = 100 satuan  waktu.  Simulasi  dilakukan  dengan  menvariasikan  nilai  beberapa  parameter yang mempengaruhi model.
1. Kasus pertama
Untuk  simulasi  pertama  diberikan  nilai  parameter  yang  menghasilkan  populasi prey terinfeksi  punah  namun  populasi predator dan prey rentan  ada  (hidup  berdampingan). Nilai parameter yang diberikan 𝑟 = 2, 𝑎 2 = 1, 𝑐 = 0,3, 𝑐 1 = 1, 𝑐 2 = 1, 𝐾 = 3, 𝐾 1 = 0,6, 𝐾 2 = 0,5, 𝛽 = 0,3 Hasil  analitik  dari  nilai  parameter  kasus  pertama  menunjukkan kepadatan populasi prey rentan 3%, prey terinfeksi 0% dan predator 0,5% . Simulasi 1 menggunakan  metode  Adams  Bashforth  Moulton  orde  empat.  Pada  simulasi  ini menunjukkan  kepadatan  populasi prey rentan  terlihat  naik  hingga 3 pada 𝑡 = 32 kemudian stabil pada 𝑡 berikutnya karena tidak adanya predasi dan laju infeksi yang kecil. Kepadatan  populasi prey terinfeksi  turun  hingga 0 pada 𝑡 = 32 kemudian  stabil  pada 𝑡 berikutnya.  Populasi prey terinfeksi  mengalami  kepunahan  karena  laju  infeksi  yang rendah  dan  adanya  predasi.  Kepadatan  populasi predator bergerak  menurun  hingga 0,5 pada 𝑡 = 34 kemudian stabil pada 𝑡 berikutnya. Grafik predator mengikuti pergerakan grafik  populasi prey terinfeksi  karena prey terinfeksi  sebagai  sumber  makanan  juga bergerak  menurun.  Populasi predator tidak  punah  karena predator mempunyai kemungkinan  untuk  tidak  hanya  bergantung  pada prey terinfeksi  saja  (dilambangkan dengan  adanya  daya  dukung  lingkungan  yaitu 𝐾 2 ).  Disini  populasi  predator  dan  prey rentan  mengikuti  titik  kesetimbangan  ( 𝐸 4 )  yaitu  dengan  syarat  stabil  asimtotik 𝛽 < 𝑐𝑐 2 𝐾 1 +𝑐 1 𝑎 2 𝐾 2 𝐾𝑐 2 𝐾 1 . Didapatkan 𝛽 = 0,3 dan 𝑐𝑐 2 𝐾 1 +𝑐 1 𝑎 2 𝐾 2 𝐾𝑐 2 𝐾 1 = 0,38 sehingga 0,3 < 0,38 .
2. Kasus kedua
Untuk simulasi kedua diberikan nilai parameter yang menghasilkan kepadatan populasi predator, prey rentan  dan prey terinfeksi sama-sama  ada  (hidup  berdampingan).  Nilai parameter  yang  diberikan 𝑟 = 2, 𝑎 2 = 1, 𝑐 = 0.3, 𝑐 1 = 1, 𝑐 2 = 1, 𝐾 = 3, 𝐾 1 = 0.6, 𝐾 2 = 0.5, 𝛽 = 1.1 Hasil  analitik  dari  nilai  parameter  kasus  pertama  menunjukkan kepadatan  populasi prey rentan  1,11%, prey terinfeksi  0,71%  dan predator 1,2% . Simulasi  2  menunjukan  ketiga  populasi  mengalami  fluktuasi  beberapa  saat  kemudian stabil pada t berikutnya dan hidup berdampingan. Kepadatan populasi prey rentan mulai stabil  pada 𝑡 = 45 dengan  kepadatan 1,11 .  Kepadatan  populasi prey terinfeksi  mulai stabil pada 𝑡 = 45 dengan kepadatan 0,71 . Kepadatan populasi Predator mulai stabil pada 𝑡 = 41 dengan kepadatan 1,2 .
Pada simulasi 2 ketiga populasi ini mengikuti syarat titik kesetimbangan dimana ketiga populasi tersebut hidup berdampingan (𝐸 * ) yaitu sebagai berikut:
𝑐 2 𝐼 * +𝐾 2 > 𝑐 1 𝐼 * (𝐼 * +𝐾 1 ) 2 ⟺ 0,8 > 0,4 ,
𝐾 1 > 𝐾 2 ⇔ 0,6 > 0,5 .
3. Kasus ketiga
Simulasi  ketiga  ini  untuk  mengetahui  pengaruh  laju  infeksi  ( 𝛽) pada  model.  Nilai parameter yang diberikan 𝑟 = 2, 𝐾 = 3, 𝑐 = 0,3, 𝑐 1 = 1, 𝑐 2 = 1, 𝐾 1 = 0,6,
𝐾 2 = 0,5, 𝑎 2 = 1, dengan variasi beberapa nilai 𝛽 yaitu 𝛽 = 0,1, 𝛽 = 1,5, 𝛽 = 8 . (a) variasi 𝛽 = 0,1 (b) variasi 𝛽 = 1,5 (c) variasi 𝛽 = 8
Gambar 1. Grafik variasi nilai 𝛽
Pada Gambar 1 dengan variasi nilai 𝛽 dapat disimpulkan bahwa semakin besar laju infeksi maka kepadatan populasi prey terinfeksi akan berangsur naik, sehingga kepadatan prey rentan  berangsur  turun.  Kepadatan  populasi predator akan  dapat  mempertahankan kepadatan populasinya karena makanannya yaitu prey terinfeksi ada/naik.
4. Kasus keempat
Simulasi  keempat  ini  bertujuan  untuk  mengetahui  pengaruh  nilai  maksimum  laju  per kapita  pengurangan 𝐼 karena 𝑦(𝑐 1 ) pada  sistem.  Nilai  parameter 𝑟 = 2, 𝐾 = 3, 𝑐 = 0,3, 𝛽 = 2, 𝑐 2 = 1, 𝐾 1 = 0,6, 𝐾 2 = 0,5 𝑎 2 = 1, 𝑆(0) = 1, 𝐼(0) = 0,5, 𝑦(0) = 1,2 dengan variasi beberapa nilai 𝑐 1 yaitu 𝑐 1 = 0,1, 𝑐 1 = 3, 𝑐 1 = 8 .
(a) variasi 𝑐 1 = 0,1 (b) variasi 𝑐 1 = 3
(c) variasi 𝑐 1 = 8 Gambar  2.  Variasi nilai 𝑐 1
Dari Gambar 2 variasi nilai 𝑐 1 dapat disimpulkan bahwa semakin besar nilai maksimum laju  per  kapita  pengurangan I karena 𝑦(𝑐 1 ) ,  kepadatan prey terinfeksi  akan  berangsur turun kemudian berfluktuasi dan berakhir punah. Kepadatan populasi prey rentan akan naik  karena prey terinfeksi  akan  berangsur  punah  dengan  adanya  predasi.  Kepadatan populasi predator akan menurun karena makanannya yaitu prey terinfeksi menurun dan punah.
5. Kasus kelima
Simulasi kelima ini bertujuan untuk mengetahui pengaruh nilai maksimum laju per kapita pengurangan 𝑦 karena 𝐼 (𝑐 2 ) pada  sistem.  Adapun  nilai  parameter 𝑟 = 2, 𝐾 = 3, 𝑐 = 0,3, 𝑐 1 = 1, 𝛽 = 2, 𝐾 1 = 0,6, 𝐾 2 = 0,5, 𝑎 2 = 1, 𝑆(0) = 1, 𝐼(0) = 0,5, 𝑦(0) = 1,2 dengan variasi beberapa nilai 𝑐 2 yaitu 𝑐 2 = 0,1, 𝑐 2 = 0,2, 𝑐 2 = 0,6 .
Gambar 3 variasi  nilai 𝑐 2 dengan nilai kurang dari satu. Pada Gambar 5 (a) kepadatan populasi predator mencapai  5%,  yang mana itu tidak dapat  diterima karena kepadatan populasi predator tidak  bisa  berkembang,  populasi predator hanya  bergantung  pada makannya dan daya dukung lingkungan saja . Ketika nilai 𝑐 2 mendekati angka satu maka kepadatan populasi predator akan mencapai 1.2 (seperti nilai awal). Semakin besar nilai 𝑐 2 (lebih dari satu) akan membuat populasi predator menurun. Kepadatan populasi prey rentan akan turun. Kepadatan populasi prey terinfeksi akan naik dan ketiga populasi bisa hidup berdampingan.
(a) variasi 𝑐 2 = 0.1 (b) variasi 𝑐 2 = 0.2 (c) variasi 𝑐 2 = 0.6
Gambar 3.  Grafik variasi nilai 𝑐 2

## 4. Kesimpulan
Kesimpulan yang dapat diambil dari tulisan ini adalah sebagai berikut:
1. Model  yang  diteliti  mempunyai  enam  titik  kesetimbangan  yaitu 𝐸 1 , 𝐸 2 , 𝐸 3 , 𝐸 4 , 𝐸 5 𝑑𝑎𝑛 𝐸 * .  Kemudian  dianalisis  kestabilannya  menghasilkan 𝐸 1 , 𝐸 2 , 𝐸 3 , 𝐸 5 tidak stabil. 𝐸 4 dan 𝐸 * adalah stabil asimtotik dengan syarat.
2. Metode Adams Basforth Moulton orde empat menghasilkan penyelesaian yang cukup akurat dan efisien.
3. Hasil interpretasi simulasi dengan variasi beberapa parameter adalah sebagai berikut:
a) Semakin besar laju infeksi (𝛽) pada prey maka populasi prey terinfeksi mengalami kenaikan  yang  menyebabkan  populasi prey rentan  menurun.  Populasi predator meningkat karena prey terinfeksi sebagai makanannya mengalami kenaikan.
b) Semakin besar nilai maksimum laju per kapita pengurangan prey terinfeksi karena predator (𝑐 1 ) maka  populasi prey terinfeksi  akan  menurun  kemudian  punah sehingga populasi prey rentan akan naik. Populasi predator akan menurun karena prey terinfeksi sebagai makanannya menurun.

## Daftar Pustaka
[1] Bacaer,  N.  (2011). A  Short  History  of  Mathematical  Population  Dynamics. New York (US): Springer-Verlag.
[2] Boyce,  W.  E.  dan  Diprima,  R.  C.  (2012). Elementary  Differensial  Equations  and Boundary Value Problems. Tenth Edition . John Willey & Sons Inc. New York.
[3] Chattopadhyay, J., Bairagi N., Chaudhuri S.. (2007). Harvesting as disease control measure  in  an  eco-epidemiological  system-atheoritical  study. Mathematical Biosciences. 217(2): 14-144.
[4] Leslie, P. H. dan Gower, J. C. (1960). The Properties of a Stochastic Model for the Predator-Prey Type of Interaction between Two Species. Biometrika , 47 (3/4): 219- 234.
[5] Mardiana,  R.  (2015).  Analisis  Solusi  Numerik  Model  Predator  Prey  Tipe  Holling dengan  Faktor  Pemanenan  pada  Prey  Menggunakan  Metode  Adams  Bashforth- Moulton. Skripsi. Jember: Universitas Jember.
[6] Putri,  P.  P.  (2013).  Analisis  Solusi  Numerik  Model  Predator-Prey  dengan  Metode Runge-Kutta Orde Empat dan Gill. Skripsi. Jember: Universitas Jember.
[7] Sahid. (2006). Pengantar Komputasi Numerik dengan MATLAB. Yogyakarta: Andi
[8] Triatmodjo, B. (1992). Metode Numerik. Yogyakarta: Beta Offset.
[9] Zhou,  Xueyoung,  Cui,  J.,  Shi,  X.,  Song,  X.  (2010).  A  Modified  Leslie-Gower predator-prey model with prey infection. J. Appl. Math. Comput. 33:471-487.