## Contemporary Mathematics and Applications
Vol. 3, No. 2, 2021, pp. 88-99.
Analisis Kestabilan Model Predator-Prey dengan Adanya Faktor Tempat Persembunyian Menggunakan Fungsi Respon Holling Tipe III
Riris Nur Patria Putri 1 , Windarto 2,* , & Cicik Alfiniyah 3
Departemen Matematika, Fakultas Sains dan Teknologi, Universitas Airlangga
*Corresponding author: windarto@fst.unair.ac.id
Abstrak. Predation  is  interaction  between  predator  and  prey,  where  predator  preys  prey.  So predators can grow, develop, and reproduce. In order for prey to avoid predators, then prey needs a refuge.  In  this  thesis,  a  predator-prey  model  with  refuge  factor  using  Holling  type  III  response function which has three populations, i.e. prey population in the refuge, prey population outside the refuge, and predator population. From the model, three equilibrium points were obtained, those are extinction of the three populations which is unstable, while extinction of predator population and coexistence are asymptotic stable under certain conditions. The numerical simulation results show that refuge have an impact the survival of the prey.
Keywords: Predator-prey model, Refuge, Holling type III, Response Function, Stability

## 1 Pendahuluan
Biosfer  merupakan  hal  yang  penting  bagi  kehidupan  makhluk  hidup,  yang  mana mencakup komponen biotik  dan abiotik.  Komponen abiotik  terdiri dari udara, daratan, air,  cahaya,  iklim,  cuaca,  suhu,  angin,  kelembapan,  kandungan  mineral,  kandungan keasaman, dan salinitas. Sedangkan komponen biotik mencakup semua makhluk hidup yang  tinggal  di  suatu  lingkungan,  seperti  manusia,  hewan,  tumbuhan,  dan mikroorganisme. Dengan adanya komponen-komponen tersebut, sangat memungkinkan terjadi  interaksi  di  dalamnya.  Salah  satunya  adalah  interaksi  antar  sesama  komponen biotik, contohnya adalah predasi [1].
Predasi  merupakan  interaksi  antara  pemangsa  ( predator )  dan  mangsa  ( prey ),  dimana predator memangsa prey .  Dalam  predasi, predator mengambil  energi  dari prey untuk tumbuh  dan  bereproduksi,  sedangkan prey merupakan  organisme  yang  dirugikan  atau dikonsumsi  oleh predator [2].  Jika  hal  ini  terjadi  terus-menerus,  maka prey akan mengalami kepunahan. Adapun cara supaya prey terhindar dari kepunahan, salah satunya dengan  mencari  tempat  persembunyian.  Tempat  persembunyian  dapat  mempengaruhi dinamika  populasi predator - prey .  Tempat  persembunyian prey memiliki  dampak  yang sangat besar terhadap koeksistensi predator dan prey dengan cara mengurangi kepunahan pada prey . Jenis ekosistem ini dapat ditemukan pada ekosistem terumbu karang [3].
Penelitian pemodelan matematika mengenai predator-prey telah banyak dilakukan. Kar [4] memodelkan dan menganalisa sistem prey - predator dengan memperhatikan adanya
tempat persembunyian mangsa dan pemanenan di kedua spesies. Dubey [5] memodelkan model predator - prey dengan  memperhatikan  adanya  area  persembunyian  yang dilindungi.  Wang  dkk  [6]  memodelkan  piramida  biomassa  terbalik  dan  tempat persembunyian  pada  ekosistem.  Berezovskaya  dkk  [7]  mengamati  peran  penyebaran mangsa  dan  tempat  persembunyian  pada  dinamika predator - prey .  Singh  dkk  [8] memodelkan  struktur  biomassa  ikan  pada  terumbu  karang  yang  hampir  murni  dan degradasi karena adanya penangkapan ikan. Wang dkk [3] meneliti dinamika predator - prey dan piramida biomassa dengan faktor tempat persembunyian.
Dari  beberapa  alasan  di  atas,  penulis  ingin  melihat  bagaimana  model predator - prey dengan adanya faktor tempat persembunyian menggunakan fungsi respon Holling tipe III dari waktu ke waktu. Selain itu, penulis ingin mengetahui bagaimana titik setimbang yang dihasilkan  dari  model predator - prey tersebut,  dan  juga  ingin  mengetahui  bagaimana kestabilan  dari  model predator - prey dengan  adanya  faktor  tempat  persembunyian menggunakan fungsi respon Holling tipe III di sekitar titik setimbangnya. Adapun batasan masalah yang penulis gunakan dalam penulisan ini yakni penulis hanya memperhatikan model predator - prey pada  terumbu  karang  yang  ditulis  oleh  Wang  dkk  [3].  Sehingga penulis tidak memperhatikan model predator - prey populasi hewan lainnnya.

## 2 Model Matematika
Model dasar predator-prey ini dirujuk dari artikel ilmiah Wang dkk. [3]. Model predator - prey ini  memiliki  tiga  kompartemen,  yaitu 𝑥 𝑟 , 𝑥 ,  dan 𝑦 dimana  secara  berurutan menunjukkan  kepadatan  populasi prey yang  berada  di  dalam  tempat  perlindungan, kepadatan populasi prey yang berada di luar tempat perlindungan, dan kepadatan populasi predator . Dalam kasus ini terdapat beberapa asumsi sebagai berikut:
a. Prey hanya  dapat  menemukan  makanan  dan  tumbuh  di  dalam  tempat perlindungan,  dan  kapasitas  pendukung  merupakan  fungsi  naik  dari  ukuran tempat persembunyian.
b. Pertumbuhan prey di  dalam  tempat  perlindungan  merupakan  pertumbuhan logistik.
c. Pada saat ketiadaan prey , predator turun secara eksponensial.
d. Predator hanya dapat memangsa prey yang berada di luar tempat perlindungan.
e. Prey dapat  migrasi  keluar masuk tempat  perlindungan, dan hanya laju imigrasi yang bergantung terhadap ukuran tempat perlindungan.
f. Laju kematian prey yang berada di luar tempat persembunyian diabaikan.
Variabel dan parameter yang akan digunakan pada model predator-prey dengan adanya faktor tempat persembunyian menggunakan fungsi respon holling tipe III akan diberikan sebagai berikut:
Tabel 1 Definisi Variabel Model Predador - Prey Variabel Keterangan Satuan 𝑥 𝑟 Kepadatan populasi prey yang berada di dalam tempat persembunyian
𝑘𝑔 𝑚 3 ⁄
𝑥
Kepadatan populasi prey yang berada di luar tempat persembunyian
𝑘𝑔 𝑚 3 ⁄
𝑦
Kepadatan populasi predator 𝑘𝑔 𝑚 3 ⁄
Tabel 2 Definisi Parameter Model Predador - Prey Parameter Keterangan Satuan 𝑎 Laju pertumbuhan per kapita dari prey yang berada di dalam tempat persembunyian 1 ℎ𝑎𝑟𝑖 ⁄ 𝐾(𝑟) Kapasitas pendukung prey di dalam tempat persembunyian yang bergantung terhadap ukuran tempat persembunyian dengan 𝐾(𝑟) = 𝐾𝑟 𝑘𝑔 𝑚 3 ⁄
𝑟 Ukuran tempat persembunyian prey - 𝛼 Laju emigrasi dari tempat persembunyian 1 ℎ𝑎𝑟𝑖 ⁄ 𝛽(𝑟) Laju imigrasi ke tempat persembunyian yang bergantung terhadap ukuran tempat persembunyian dengan 𝛽(𝑟) = 𝛽𝑟 2 𝜉 2 +𝑟 2 1 ℎ𝑎𝑟𝑖 ⁄ 𝛾 Laju maksimum predasi 𝑚 6 𝑘𝑔 2 ⁄ ℎ𝑎𝑟𝑖 𝜃
Tingkat kejenuhan predator 𝑚 6 𝑘𝑔 2 ⁄ 𝑐 Nilai konversi prey yang diubah menjadi predator - 𝑑 Laju kematian per kapita predator 1 ℎ𝑎𝑟𝑖 ⁄
Pada  model predator-prey dengan  adanya  faktor  tempat  persembunyian  menggunakan fungsi respon Holling tipe III diasumsikan nilai variabel 𝑥 𝑟 , 𝑥, 𝑦  ≥ 0 dan nilai parameter bernilai positif.
Berikut adalah model predator-prey dengan adanya faktor tempat persembunyian yang digunakan pada penulisan ini:
𝑑𝑥 𝑟 𝑑𝑡 = 𝑎𝑥 𝑟 [1 − 𝑥 𝑟 𝐾𝑟 ] − 𝛼𝑥 𝑟 + 𝛽𝑟 2 𝜉 2 +𝑟 2 𝑥 (1)
𝑑𝑥 𝑑𝑡 = 𝛼𝑥 𝑟 − 𝛽𝑟 2 𝜉 2 +𝑟 2 𝑥 − 𝛾𝑥 2 1+𝜃𝑥 2 𝑦 (2)
𝑑𝑦 𝑑𝑡 = 𝑐 𝛾𝑥 2 1+𝜃𝑥 2 𝑦 − 𝑑𝑦 (3)

## 3 Analisis Kestabilan Model
Pada  bagian  ini  akan  dibahas  megenai  analisis  kestabilan  lokal  dari  titik  setimbang kepunahan  ketiga  populasi,  titik  setimbang  kepunahan predator ,  serta  titik  setimbang koeksistensi.

## 3.1. Titik Setimbang
Berdasarkan Abell dan Braselton [9], titik setimbang merupakan titik pada saat kondisi ketika perubahan jumlah populasi tertentu sepanjang satuan waktu adalah nol. Keadaan setimbang pada model predator-prey ini terpenuhi saat 𝑑𝑥 𝑟 𝑑𝑡 = 𝑑𝑥 𝑑𝑡 = 𝑑𝑦 𝑑𝑡 = 0 .
Titik  setimbang  kepunahan  ketiga  populasi merupakan  titik  setimbang  pada  saat kondisi populasi prey yang berada di dalam tempat persembunyian, populasi prey yang berada  di  luar  tempat  persembunyian,  dan  populasi predator sepanjang  satuan  waktu mengalami  kepunahan.  Kondisi  ini  artinya 𝑥 𝑟 = 0 , 𝑥 = 0 , 𝑦 = 0 .  Oleh  karena  itu, diperoleh titik setimbang 𝐸 0 = (𝑥 𝑟0 , 𝑥 0 , 𝑦 0 ) = (0, 0, 0) .
Titik  setimbang  kepunahan predator merupakan  titik  setimbang  pada  saat  kondisi populasi predator mengalami  kepunahan,  sehingga  tidak  adanya predator yang memangsa prey yang  berada  di  luar  tempat  persembunyian.  Kondisi  ini  terjadi  ketika 𝑥 𝑟 ≠ 0 , 𝑥 ≠ 0 , 𝑦 = 0 .  Sehingga  dapat  dimisalkan  titik  setimbang  kepunahan predator adalah 𝐸 1 = (𝑥 𝑟1 , 𝑥 1 , 𝑦 1 ) = (𝑥 𝑟1 , 𝑥 1 , 0) .  Dari  persamaan  (1),  (2),  (3)  diperoleh 𝑑𝑥 𝑟 𝑑𝑡 = 𝑑𝑥 𝑑𝑡 = 𝑑𝑦 𝑑𝑡 = 0 .  Sehingga  diperoleh  titik  setimbang  kepunahan predator 𝐸 1 = (𝑥 𝑟1 , 𝑥 1 , 𝑦 1 ) = (𝐾𝑟, 𝛼𝐾𝑟(𝜉 2 +𝑟 2 ) 𝛽𝑟 2 , 0) .
Titik  setimbang  koeksistensi merupakan  titik  setimbang  ketika  kondisi  populasi prey yang berada di dalam tempat persembunyian, populasi prey yang berada di luar tempat persembunyian, dan populasi predator dapat hidup berdampingan dalam suatu ekosistem. Kondisini ini terjadi ketika 𝑥 𝑟 ≠ 0 , 𝑥 ≠ 0 , 𝑦 ≠ 0 . Dari persamaan (1), (2), (3) didapatkan 𝑑𝑥 𝑟 𝑑𝑡 = 𝑑𝑥 𝑑𝑡 = 𝑑𝑦 𝑑𝑡 = 0 sehingga  titik  setimbang  koeksistensi 𝐸 2 = (𝑥 𝑟2 , 𝑥 2 , 𝑦 2 ) dengan didapat sebagai berikut:
𝑥 2 = √ 𝑑 𝑐𝛾−𝑑𝜃 (4)
𝑥 𝑟2 = 𝐾𝑟 [ (𝑎−𝛼)+√(𝑎−𝛼) 2 +4 𝑎𝛽𝑟
𝐾(𝜉2+𝑟2) 𝑥 2 2𝑎 ] (5)
𝑦 2 = 𝑐 𝑑 ( 𝛼𝐾𝑟 2𝑎 [(𝑎 − 𝛼) + √(𝑎 − 𝛼) 2 + 4 𝑎𝛽𝑟 𝐾(𝜉 2 +𝑟 2 ) 𝑥 2 ] − 𝛽𝑟 2 𝜉 2 +𝑟 2 𝑥 2 ) (6)
Titik setimbang koeksistensi ada
jika 𝑐𝛾 − 𝑑𝜃 > 0 , 𝑎 − 𝛼 > 0 , dan
𝛼𝐾𝑟 [ (𝑎−𝛼)+√(𝑎−𝛼) 2 +4 𝑎𝛽𝑟
𝐾(𝜉2+𝑟2) 𝑥 2
2𝑎 ] − 𝛽𝑟 2 𝜉 2 +𝑟 2 𝑥 2 > 0 .

## 3.2. Analisis Kestabilan di Sekitar Titik Setimbang
Pada  bagian  ini  akan  dilakukan  analisis  kestabilan  lokal  dari  masing-masing  titik setimbang.  Analisis  ini  dilakukan  untuk  mengetahui  dinamika  perilaku  sistem  pada model. Analisis kestabilan lokal didekati dengan menggunakan nilai eigen (𝜆) , sehingga kestabilan  asimtotis  lokal  dari  model predator - prey dengan  adanya  faktor  tempat persembunyian menggunakan fungsi respon Holling tipe III di titik setimbang 𝐸 0 , 𝐸 1 , dan 𝐸 2 dilakukan dengan menggunakan matriks  Jacobian. Dari persamaan (1), (2), dan (3) didapatkan Matriks Jacobian sebagai berikut:
𝐽 =
[ 𝑎 − 𝛼 − 2𝑎𝑥 𝑟 𝐾𝑟 𝛽𝑟 2 𝜉 2 +𝑟 2 0
𝛼 − 𝛽𝑟 2 𝜉 2 +𝑟 2 − 2𝛾𝑥𝑦[(1+𝜃𝑥 2 )−𝜃𝑥 2 𝑦] (1+𝜃𝑥 2 ) 2 − 𝛾𝑥 2 1+𝜃𝑥 2 0 2𝑐𝛾𝑥𝑦[(1+𝜃𝑥 2 )−𝜃𝑥 2 𝑦] (1+𝜃𝑥 2 ) 2 𝑐𝛾𝑥 2 1+𝜃𝑥 2 − 𝑑]
Selanjutnya,  akan dianalisis  kestabilan dari titik setimbang kepunahan ketiga populasi, kepunahan predator , dan koeksistensi.
a. Kestabilan lokal titik setimbang kepunahan ketiga populasi
Matriks Jacobian untuk titik setimbang kepunahan ketiga populasi 𝐸 0 = (𝑥 𝑟0 , 𝑥 0 , 𝑦 0 ) = (0, 0, 0) adalah sebagai berikut.
𝐽(𝐸 0 ) =
[ 𝑎 − 𝛼 𝛽𝑟 2 𝜉 2 +𝑟 2 0 𝛼 − 𝛽𝑟 2 𝜉 2 +𝑟 2 0 0 0 −𝑑]
(7)
Kemudian  matriks 𝐽(𝐸 0 ) di  atas  dapat  dibentuk  persamaaan  karakteristik  dengan menggunakan 𝑑𝑒𝑡(𝜆𝐼 − 𝐽(𝐸 0 )) = 0 sebagai berikut.
⇔ (𝜆 + 𝑑) [(𝜆 − 𝑎 + 𝛼) (𝜆 + 𝛽𝑟 2 𝜉 2 +𝑟 2 ) − 𝛼𝛽𝑟 2 𝜉 2 +𝑟 2 ] = 0 (8)
Berdasarkan persamaan (8), diperoleh nilai eigen 𝜆 1 = −𝑑 dan akar-akar dari persamaan berikut.
𝜆 2 + 𝑏 1 𝜆 + 𝑏 2 = 0 (9)
𝑏 1 = 𝛼 − 𝑎 + 𝛽𝑟 2 𝜉 2 +𝑟 2 𝑏 2 = − 𝑎𝛽𝑟 2 𝜉 2 +𝑟 2
Titik setimbang kepunahan ketiga populasi tersebut akan stabil asimtotis jika dan hanya jika persamaan (9) mempunyai akar-akar bagian real negatif. Dari penjelasan sebelumnya jelas  bahwa 𝜆 1 bernilai  negatif  sebab  semua  parameter  pada  model  berniali  positif. Selanjutnya, akan ditentukan syarat agar persamaan (9) memiliki akar-akar yang bagian realnya negatif, yaitu dengan menggunakan kriteria Routh Hurwitz.
Persamaan karakteristik (9) akan memiliki akar-akar yang bagian realnya negatif jika nilai dari 𝑏 1 , 𝑏 2 > 0 .
i. Koefisien 𝑏 1 akan bernilai lebih dari nol apabila memenuhi 𝛼 + 𝛽𝑟 2 𝜉 2 +𝑟 2 > 𝑎 . ii. Karena semua nilai parameter bernilai positif, maka dapat dipastikan bahwa 𝑏 2 < 0 .
Berdasarkan uraian di atas diperoleh bahwa titik setimbang kepunahan ketiga populasi 𝐸 0 tidak stabil.

## b. Kestabilan lokal titik setimbang kepunahan predator
Matriks  Jacobian  untuk  titik  setimbang  kepunahan predator 𝐸 1 = (𝑥 𝑟1 , 𝑥 1 , 0) adalah sebagai berikut.
𝐽(𝐸 1 ) = [ −𝑎 − 𝛼 𝛽𝑟 2 𝜉 2 +𝑟 2 0 𝛼 − 𝛽𝑟 2 𝜉 2 +𝑟 2 − 𝛾(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 (𝛽𝑟 2 ) 2 +𝜃(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2
0 0
𝑐𝛾(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 (𝛽𝑟 2 ) 2 +𝜃(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 − 𝑑 ] (10)
Kemudian  matriks 𝐽(𝐸 1 ) di  atas  dapat  dibentuk  persamaaan  karakteristik  dengan menggunakan 𝑑𝑒𝑡(𝜆𝐼 − 𝐽(𝐸 1 )) = 0 sebagai berikut.
⇔ [𝜆 − 𝑐𝛾(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 (𝛽𝑟 2 ) 2 +𝜃(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 + 𝑑] [(𝜆 + 𝑎 + 𝛼) (𝜆 + 𝛽𝑟 2 𝜉 2 +𝑟 2 ) − 𝛼𝛽𝑟 2 𝜉 2 +𝑟 2 ] = 0 (11)
Berdasarkan persamaan (11), diperoleh nilai eigen 𝜆 1 = 𝑐𝛾(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 (𝛽𝑟 2 ) 2 +𝜃(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 − 𝑑 dan akar-akar dari persamaan berikut.
𝜆 2 + 𝑏 1 𝜆 + 𝑏 2 = 0 (12)
𝑏 1 = 𝛼 + 𝑎 + 𝛽𝑟 2 𝜉 2 +𝑟 2
𝑏 2 = 𝑎𝛽𝑟 2 𝜉 2 +𝑟 2
Titik  setimbang  kepunahan predator tersebut  akan  stabil  asimtotis  jika  dan  hanya  jika persamaan (12) mempunyai akar-akar bagian real negatif. Karena semua nilai parameter bernilai positif, maka dari penjelasan sebelumnya nilai 𝜆 1 akan bernilai negatif apabila 𝑐𝛾(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2
(𝛽𝑟 2 ) 2 +𝜃(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 < 𝑑 .  Selanjutnya,  akan  ditentukan  syarat  agar  persamaan  (12) memiliki  akar-akar  yang  bagian  realnya  negatif,  yaitu  dengan  menggunakan  kriteria Routh Hurwitz.
Persamaan karakteristik  (12) akan memiliki akar-akar  yang bagian realnya negatif jika nilai dari 𝑏 1 , 𝑏 2 > 0 . i. Karena semua nilai parameter bernilai positif, maka dapat dipastikan bahwa 𝑏 1 >
0 .
ii. Karena semua nilai parameter bernilai positif, maka dapat dipastikan bahwa 𝑏 2 > 0 .
Berdasarkan uraian di atas diperoleh bahwa titik setimbang kepunahan predator 𝐸 1 stabil asimtotis jika 𝑐𝛾(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 (𝛽𝑟 2 ) 2 +𝜃(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 < 𝑑 .

## c. Kestabilan lokal titik setimbang koeksistensi
Matriks  Jacobian  untuk  titik  setimbang  koeksistensi 𝐸 2 = (𝑥 𝑟2 , 𝑥 2 , 𝑦 2 ) adalah  sebagai berikut.
𝐽(𝐸 2 ) =
[ 𝑎 − 𝛼 − 2𝑎𝑥 𝑟 𝐾𝑟 𝛽𝑟 2 𝜉 2 +𝑟 2 0 𝛼 − 𝛽𝑟 2 𝜉 2 +𝑟 2 − 2𝛾𝑥𝑦[(1+𝜃𝑥 2 )−𝜃𝑥 2 𝑦] (1+𝜃𝑥 2 ) 2 − 𝛾𝑥 2 1+𝜃𝑥 2
0 2𝑐𝛾𝑥𝑦[(1+𝜃𝑥 2 )−𝜃𝑥 2 𝑦] (1+𝜃𝑥 2 ) 2 𝑐𝛾𝑥 2 1+𝜃𝑥 2 − 𝑑]
(13)
Dengan menggunakan kondisi kesetimbangan koeksistensi di dapat matriks 𝐽(𝐸 2 ) sebagai berikut:
𝐽(𝐸 2 ) =
[ 𝑎 − 𝛼 − 2𝑎𝑥 𝑟 𝐾𝑟 𝛽𝑟 2 𝜉 2 +𝑟 2 0 𝛼 − 𝛽𝑟 2 𝜉 2 +𝑟 2 − 2𝛾𝑥𝑦[(1+𝜃𝑥 2 )−𝜃𝑥 2 𝑦] (1+𝜃𝑥 2 ) 2 − 𝑑 𝑐
0 2𝑐𝛾𝑥𝑦[(1+𝜃𝑥 2 )−𝜃𝑥 2 𝑦] (1+𝜃𝑥 2 ) 2 0 ] (14)
Kemudian  matriks 𝐽(𝐸 2 ) di  atas  dapat  dibentuk  persamaaan  karakteristik  dengan menggunakan 𝑑𝑒𝑡(𝜆𝐼 − 𝐽(𝐸 2 )) = 0 sebagai berikut: 𝜆 3 + 𝑏 1 𝜆 2 + 𝑏 2 𝜆 + 𝑏 3 = 0 (15)
𝑏 1 = 1 𝐾𝑟𝑥 𝑟 𝑥(𝜉 2 + 𝑟 2 )(1 + 𝜃𝑥 2 ) 2 (𝐾𝑎𝑟 3 𝜃 2 𝑥 𝑟 2 𝑥 4 + 𝐾𝛽𝑟 3 𝜃 2 𝑥 6 − 𝐾𝑟 3 𝛾𝜃𝑥 𝑟 𝑥 4 𝑦 + 2𝐾𝑎𝑟 3 𝜃𝑥 𝑟 2 𝑥 2 + 2𝐾𝑟 3 𝛽𝜃𝑥 4 + 𝐾𝑟 3 𝛾𝑥 𝑟 𝑥 2 𝑦 + 𝐾𝑎𝑟 3 𝑥 𝑟 2 + 𝐾𝛽𝑟 3 𝑥 2 + 𝑎𝑟 2 𝜃 2 𝑥 𝑟 2 𝑥 5 + 2𝑎𝑟 2 𝜃𝑥 𝑟 2 𝑥 3 + 𝑎𝑟 2 𝑥 𝑟 2 𝑥 + 2𝐾𝑎𝑟𝜃𝜉 2 𝑥 𝑟 2 𝑥 2 + 𝐾𝑎𝑟𝜃 2 𝜉 2 𝑥 𝑟 2 𝑥 4 − 𝐾𝑟𝛾𝜃𝜉 2 𝑥 𝑟 𝑥 4 𝑦 + 𝐾𝑟𝛾𝜉 2 𝑥 𝑟 𝑥 2 𝑦 + 𝐾𝑎𝑟𝜉 2 𝑥 𝑟 2 + 𝑎𝜃 2 𝜉 2 𝑥 𝑟 2 𝑥 5 + 2𝑎𝜃𝜉 2 𝑥 𝑟 2 𝑥 3 + 𝑎𝜉 2 𝑥 𝑟 2 𝑥) 𝑏 2 = 1 𝐾𝑟𝑥 𝑟 𝑥(𝜉 2 + 𝑟 2 )(1 + 𝜃𝑥 2 ) 2 (𝐾𝑎𝑟 3 𝛽𝜃 2 𝑥 𝑟 𝑥 5 − 𝐾𝑟 3 𝛼𝛽𝜃 2 𝑥 𝑟 𝑥 5 − 𝐾𝑟 3 𝛽𝛾𝜃𝑥 5 𝑦 + 2𝐾𝑎𝑟 3 𝛽𝜃𝑥 𝑟 𝑥 3 − 2𝐾𝑟 3 𝛼𝛽𝜃𝑥 𝑟 𝑥 3 + 𝐾𝑟 3 𝛽𝛾𝑥 3 𝑦 + 2𝐾𝑑𝑟 3 𝛾𝑥 𝑟 𝑥 2 𝑦 + 𝐾𝑎𝑟 3 𝛽𝑥 𝑟 𝑥 − 𝐾𝑟 3 𝛼𝛽𝑥 𝑟 𝑥 + 𝑎 2 𝑟 2 𝜃 2 𝑥 𝑟 3 𝑥 4 − 𝑎𝑟 2 𝛾𝜃𝑥 𝑟 2 𝑥 4 𝑦 + 2𝑎 2 𝑟 2 𝜃𝑥 𝑟 3 𝑥 2 + 𝑎𝑟 2 𝛾𝑥 𝑟 2 𝑥 2 𝑦 + 𝑎 2 𝑟 2 𝑥 𝑟 3 + 2𝐾𝑑𝑟𝛾𝜉 2 𝑥 𝑟 𝑥 2 𝑦 + 𝑎 2 𝜃 2 𝜉 2 𝑥 𝑟 3 𝑥 4 − 𝑎𝛾𝜃𝜉 2 𝑥 𝑟 2 𝑥 4 𝑦 + 2𝑎 2 𝜃𝜉 2 𝑥 𝑟 3 𝑥 2 + 𝑎𝛾𝜉 2 𝑥 𝑟 2 𝑥 2 𝑦 + 𝑎 2 𝜉 2 𝑥 𝑟 3 ) 𝑏 3 = 1 𝐾𝑟𝑥 𝑟 (𝜉 2 + 𝑟 2 )(1 + 𝜃𝑥 2 ) 2 (2𝑑𝛾𝑥𝑦(𝐾𝑟 3 𝛽𝑥 + 𝑎𝑟 2 𝑥 𝑟 2 + 𝑎𝜉 2 𝑥 𝑟 2 ))
Titik  setimbang  koeksistensi  akan  stabil  asimtotis  jika  dan  hanya  jika  persamaan  (15) mempunyai  akar-akar  bagian  real  negatif.  Selanjutnya,  akan  ditentukan  syarat  agar persamaan memiliki akar-akar yang bagian realnya negatif, yaitu dengan menggunakan kriteria Routh-Hurwitz.
Persamaan karakteristik (15) akan memiliki akar-akar yang bagian real negatif jika nilai dari 𝑏 1 , 𝑏 2 , 𝑏 3 , (𝑏 1 𝑏 2 − 𝑏 3 ) > 0 . i. Agar memenuhi 𝑏 1 > 0 , maka harus memenuhi syarat adalah sebagai berikut: 𝐾𝑎𝑟 3 𝜃 2 𝑥 𝑟 2 𝑥 4 + 𝐾𝛽𝑟 3 𝜃 2 𝑥 6 + 2𝐾𝑎𝑟 3 𝜃𝑥 𝑟 2 𝑥 2 + 2𝐾𝑟 3 𝛽𝜃𝑥 4 + 𝐾𝑟 3 𝛾𝑥 𝑟 𝑥 2 𝑦 + 𝐾𝑎𝑟 3 𝑥 𝑟 2 + 𝐾𝛽𝑟 3 𝑥 2 + 𝑎𝑟 2 𝜃 2 𝑥 𝑟 2 𝑥 5 + 2𝑎𝑟 2 𝜃𝑥 𝑟 2 𝑥 3 + 𝑎𝑟 2 𝑥 𝑟
2 𝑥 + 2𝐾𝑎𝑟𝜃𝜉 2 𝑥 𝑟
2 𝑥 2 + 𝐾𝑎𝑟𝜃 2 𝜉 2 𝑥 𝑟 2 𝑥 4 + 𝐾𝑟𝛾𝜉 2 𝑥 𝑟 𝑥 2 𝑦 + 𝐾𝑎𝑟𝜉 2 𝑥 𝑟
2 + 𝑎𝜃 2 𝜉 2 𝑥 𝑟 2 𝑥 5 + 2𝑎𝜃𝜉 2 𝑥 𝑟 2 𝑥 3 + 𝑎𝜉 2 𝑥 𝑟 2 𝑥 > 𝐾𝑟 3 𝛾𝜃𝑥 𝑟 𝑥 4 𝑦 + 𝐾𝑟𝛾𝜃𝜉 2 𝑥 𝑟 𝑥 4 𝑦 ii. Syarat agar 𝑏 2 > 0 adalah sebagai berikut:
𝐾𝑎𝑟 3 𝛽𝜃 2 𝑥 𝑟 𝑥 5 + 2𝐾𝑎𝑟 3 𝛽𝜃𝑥 𝑟 𝑥 3 + 𝐾𝑟 3 𝛽𝛾𝑥 3 𝑦 + 2𝐾𝑑𝑟 3 𝛾𝑥 𝑟 𝑥 2 𝑦 + 𝐾𝑎𝑟 3 𝛽𝑥 𝑟 𝑥 + 𝑎 2 𝑟 2 𝜃 2 𝑥 𝑟 3 𝑥 4 + 2𝑎 2 𝑟 2 𝜃𝑥 𝑟 3 𝑥 2 + 𝑎𝑟 2 𝛾𝑥 𝑟 2 𝑥 2 𝑦 + 𝑎 2 𝑟 2 𝑥 𝑟 3 +
2𝐾𝑑𝑟𝛾𝜉 2 𝑥 𝑟 𝑥 2 𝑦 + 𝑎 2 𝜃 2 𝜉 2 𝑥 𝑟 3 𝑥 4 + 2𝑎 2 𝜃𝜉 2 𝑥 𝑟 3 𝑥 2 + 𝑎𝛾𝜉 2 𝑥 𝑟 2 𝑥 2 𝑦 + 𝑎 2 𝜉 2 𝑥 𝑟 3 > 𝐾𝑟 3 𝛼𝛽𝜃 2 𝑥 𝑟 𝑥 5 + 𝐾𝑟 3 𝛽𝛾𝜃𝑥 5 𝑦 + 2𝐾𝑟 3 𝛼𝛽𝜃𝑥 𝑟 𝑥 3 + 𝐾𝑟 3 𝛼𝛽𝑥 𝑟 𝑥 + 𝑎𝑟 2 𝛾𝜃𝑥 𝑟
2 𝑥 4 𝑦 + 𝑎𝛾𝜃𝜉 2 𝑥 𝑟 2 𝑥 4 𝑦 iii. Karena 𝑥 𝑟 , 𝑥, 𝑦 > 0 dan nilai parameter bernilai positif, maka 𝑏 3 > 0.
Dalam menentukan syarat (𝑏 1 𝑏 2 − 𝑏 3 ) > 0 untuk menentukan kestabilan titik setimbang koeksistensi  sulit  untuk  dilakukan  secara  analitik,  maka  dilakukan  secara  numerik
menggunakan  software  MATLAB.  Simulasi  ini  dilakukan  dengan  memberi  nilai parameter  dan  empat  nilai  awal  untuk 𝑥 𝑟 (0) , 𝑥(0) ,  dan 𝑦(0) yang  diketahui.  Berikut adalah tabel nilai awal dan nilai parameter yang digunakan. Simulasi ini hanya dilakukan pada waktu 𝑡 = 0 sampai 𝑡 = 50000 hari.
Tabel 3 Nilai Awal untuk Bidang Fase Koeksistensi
Nilai Awal 𝑥 𝑟 (0) 𝑥(0) 𝑦(0) Warna 1 0,20 0,15 0,10 Hijau 2 0,30 0,25 0,20 Biru 3 0,40 0,35 0,30 Merah Tabel 4 Nilai Parameter Model Predator - Prey Saat Kondisi Koeksistensi Parameter Nilai Parameter Sumber 𝑎 0,6 Asumsi 𝐾 0,8 Wang dkk (2018) 𝑟 1 Asumsi 𝛼 0,035 Wang dkk (2018) 𝛽 0,0119 Wang dkk (2018) 𝛾 0,112 Asumsi 𝑐 0,01 Asumsi 𝑑 0,0007 Asumsi 𝜉 1 Asumsi 𝜃 0,1 Asumsi
Berdasarkan nilai parameter pada Tabel 3, berikut ini adaah simulasi bidang fase antara kepadatan  populasi prey yang  berada  di  luar  tempat  persembunyian  dan  kepadatan populasi predator selama rentang waktu 50000 hari.
Gambar 1 Grafik Bidang Fase 𝑥(𝑡) dan 𝑦(𝑡) pada Titik Setimbang 𝐸 2
Pada Gambar 1 adalah grafik bidang fase kepadatan populasi prey yang berada di dalam tempat  persembunyian  dan  kepadatan  populasi predator pada  model predator - prey dengan adanya faktor tempat persembunyian menggunakan fungsi respon Holling tipe III dengan syarat 𝐸 2 . Berdasarkan gambar tersebut dengan ketiga nilai awal yang digunakan menunjukkan  bahwa  semua  grafiknya  cenderung  menuju  titik  setimbang  koeksistensi, yaitu  titik (0,8165 ;  0,3115) dengan  kata  lain  grafik  tersebut  konvergen.  Dengan demikian model tersebut stabil asimtotis pada titik setimbang koeksistensi.

## 3.3. Simulasi Numerik
Pada bagian ini dibahas mengenai simulasi numerik dan interpretasi model predator - prey dengan  adanya  faktor  tempat  persembunyian  menggunakan  fungsi  respon  Holling  tipe III. Simulasi ini bertujuan untuk mengetahui dinamika populasi antara populasi prey yang berada  di  dalam  tempat  persembunyian,  populasi prey yang  berada  di  luar  tempat persembunyian, dan populasi predator pada selang waktu tertentu.
Simulasi  titik  setimbang  kepunahan predator merupakan  simulasi  numerik  model predator - prey dengan adanya faktor tempat persembunyian menggunakan fungsi respon Holling tipe III yang memenuhi kondisi kestabilan titik setimbang kepunahan predator (𝐸 1 ) . Nilai parameter yang digunakan untuk simulasi disajikan pada Tabel 5 dengan nilai awal 𝑥 𝑟 (0) = 0,30 , 𝑥(0) = 0,20 , dan 𝑦(0) = 0,50 selama 500 hari.
Tabel 5 Nilai Parameter Model Predator - Prey Saat Kondisi Kepunahan Predator Parameter Nilai Parameter Sumber 𝑎 0,6 Asumsi 𝐾 0,8 Wang dkk (2018) 𝑟 1 Asumsi 𝛼 0,035 Wang dkk (2018) 𝛽 0,0119 Wang dkk (2018) 𝛾 0,0112 Wang dkk (2018) 𝑐 0,04 Wang dkk (2018) 𝑑 0,007 Asumsi 𝜉 1 Asumsi 𝜃 0,1 Asumsi
Berikut hasil simulasi yang menggambarkan dinamika populasi antara populasi prey yang berada  di  dalam  tempat  persembunyian,  populasi prey yang  berada  di  luar  tempat persembunyian, dan populasi predator :
Gambar 2 Dinamika Predator - Prey Saat Kondisi Kepunahan Predator
Simulasi titik setimbang koeksistensi merupakan simlasi numerik model predator - prey dengan adanya faktor tempat persembunyian menggunakan fungsi respon Holling tipe III yang  memenuhi  kondisi  kestabilan  titik  setimbang  koeksistensi (𝐸 2 ) .  Nilai  parameter yang digunakan untuk simulasi disajikan pada Tabel 4 dengan nilai awal 𝑥 𝑟 (0) = 0,40 , 𝑥(0) = 0,50 , dan 𝑦(0) = 0,20 selama 5000 hari. Berikut hasil simulasi yang menggambarkan dinamika populasi antara populasi prey yang berada  di  dalam  tempat  persembunyian,  populasi prey yang  berada  di  luar  tempat persembunyian, dan populasi predator :
Gambar 3 Grafik Dinamika Predator - Prey Saat Kondisi Koeksistensi

## 4 Kesimpulan
Berdasarkan hasil yang telah diperoleh pada pembahasan, dapat disimpulkan bahwa titik setimbang kepunahan bersifat tidak stabil, sedangkan titik setimbang kepunahan predator dengan syarat 𝑐𝛾(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 (𝛽𝑟 2 ) 2 +𝜃(𝛼𝐾𝑟(𝜉 2 +𝑟 2 )) 2 < 𝑑 dan titik setimbang koeksistensi dengan syarat koeksistensi 𝑐𝛾 − 𝑑𝜃 > 0 , 𝑎 − 𝛼 > 0 , 𝛼𝐾𝑟 [ (𝑎−𝛼)+√(𝑎−𝛼) 2 +4 𝑎𝛽𝑟
𝐾(𝜉2+𝑟2) 𝑥 2 2𝑎 ] − 𝛽𝑟 2 𝜉 2 +𝑟 2 𝑥 2 > 0 bersifat  stabil  asimtotis.  Hasil  simulasi  numerik  pada  saat  kesetimbangan  kepunahan predator , kepadatan populasi prey yang berada di luar tempat persembunyian meningkat. Sedangkan pada saat kesetimbangan koeksistensi, kepadatan populasi prey yang berada di  luar  tempat  persembunyian  meningkat  hingga  pada  waktu  tertentu,  lalu  menurun hingga pada waktu tertentu dan selanjutnya cenderung konstan.

## Referensi
[1] Kistinnah,  I.,  dan  Lestari,  E.S.,  2006. Biologi  Makhluk  Hidup  dan
Lingkungannya . Jakarta: Pusat Perbukuan Departemen Pendidikan Nasional
[2] Stevens, A.N.P., https://www.nature.com/scitable , 3 Oktober 2018
[3] Wang, H., Thanarajah, S., dan Gaudreau, P., 2018, Refuge-Mediated Predator- Prey Dynamics and Biomass Pyramids, Mathematical Biosciences , 298, 29-45
[4] Kar, T. K., 2006, Modelling and Analysis of A Harvested Prey–Predator System Incorporating  A  Prey  Refuge, Journal  of  Computational  and  Applied Mathematics , 185, 19-33
[5] Dubey,  B.,  2007,  A  Prey-Predator  Model  with  a  Reserved  Area, Nonlinear Analysis: Modelling and Control , 12, 279-494
[6] Wang,  H.,  Morrison,  W.,  Singh,  A.,  dan  Weiss,  H.,  2009,  Modeling  Inverted Biomass  Pyramids  and  Refuges  in  Ecosystems, Ecological  Modelling ,  220, 1376-1382
[7] Berezovskaya,  F.  S.,  Song,  B.,  dan  Castillo-Chavez,  C.,  2010,  Role  of  Prey Dispersal  and Refuges On Predator-Prey Dynamics, SIAM J. Appl. Math. , 70, 1821-1839
[8] Singh,  A.,  Wang,  H.,  Morrision,  W.M.,  dan  Weiss,  H.,  2012,  Modeling  Fish Biomass  Structure  at  Near  Pristine  Coral  Reefs  and  Degradation  by  Fishing, Journal of Biological Systems , 20, 21-36
[9] Abell,  M.  dan  Braselton,  J.,2004, Differential  Equations  with  Mathematica , Third Edition, USA: Elsevier Academic Press