## EULER: Jurnal Ilmiah Matematika, Sains dan Teknologi

## ANALISIS KESTABILAN MODEL INTERAKSI PREDATOR-PREY

## DENGAN FUNGSI RESPON MONOD-HALDANE DAN PERILAKU ANTI PEMANGSA
Muhammad Bachtiar Gaib 1* , Wahdania At. Ja ’ a 2
1,2,3 Jurusan Matematika, Universitas Negeri Gorontalo, Bone Bolango 96119, Indonesia
*Penulis Korespondensi. Email: m.tiargaib@gmail.com

## Abstrak
Artikel  ini  mengkaji  model  mangsa-pemangsa  yang  saling  berkompetisi,  dengan  menggunakan fungsi respon Monod-Haldane dan perilaku anti pemangsa. Artikel ini membahas penentuan titik tetap, analisis kestabilan titik tetap dan simulasi numerik. Diperoleh tiga titik tetap yaitu T 1 , T 2 dan T 3 dimana kestabilan titik tetap T 1 selalu bersifat sadel, kestabilan titik tetap T 2 dan T 3 akan stabil jika memenuhi syarat  parameter yang telah ditentukan. Terdapat dua kasus dalam titik tetap T 3 , dimana kasus pertama bersifat simpul stabil dan kasus kedua bersifat spiral stabil.
Kata Kunci: Model Predator-Prey; Fungsi Respon; Monod-Haldane; Anti Pemangsa; Titik Tetap

## Abstract
This  article  examines  a  competing  prey-predator  model  using  the  Monod-Haldane  response function  and  anti-predator  behavior.  This  article  discusses  equilibrium  point  determination, equilibrium point stability analysis, and numerical simulation. Obtained three equilibrium points, namely 𝑇 1 , 𝑇 2 ,  and 𝑇 3 ,  where  the  equilibrium-point 𝑇 1 is  always  saddle,  the  stability  of  the equilibrium points 𝑇 2 and 𝑇 3 will be stable if it meets the predetermined parameter requirements. There are two cases in the equilibrium point 𝑇 3 , where the first case is vertically stable and the second case is spiral stable.
Keywords: Predator-Prey  Model;  Response  Function;  Monod-Haldane;  Anti  Predators; Equilibrium Point

## 1. Pendahuluan
Dalam kehidupan setiap makhluk hidup tidak dapat terlepas dengan yang namanya interaksi. Interaksi merupakan suatu jenis tindakan yang terjadi ketika dua atau lebih objek mempengaruhi atau memiliki efek satu sama lain. Salah satu interaksi yang terjadi adalah interaksi antara mangsa dan pemangsa, yang sering disebut dengan interaksi prey-predator . Pada interaksi predator (pemangsa) mengkonsumsi prey (mangsa) agar dapat bertahan hidup dan fungsi predator terhadap prey adalah sebagai  pengendali  populasi prey [1] .  Dinamika  populasi  adalah  cabang  ilmu  kehidupan  yang mempelajari  perubahan  jangka  pendek  dan  jangka  panjang  dalam  komposisi  ukuran  dan  usia populasi, dan proses biologis serta lingkungan yang memengaruhi perubahan tersebut. Populasi yang dinamis berkaitan dengan cara populasi dipengaruhi oleh angka kelahiran dan kematian, dan oleh imigrasi  dan  emigrasi,  dan  mempelajari  topik-topik  seperti  populasi  yang  menua  atau  penurunan populasi [2] . Model matematika merupakan representasi matematika yang dihasilkan dari pemodelan Matematika.  Pemodelan  Matematika  merupakan  suatu  proses  merepresentasikan  dan  menjelaskan permasalahan pada dunia nyata ke dalam pernyataan matematis [3] .
Beberapa  penelitian  telah  dilakukan  mengenai  model  matematika  khususnya  pada  interaksi predator-prey . Dalam [4] telah dikaji sistem predator-prey model leslie-gower dengan pemanenan
secara konstan terhadap predator . Tujuan dari penelitian tersebut adalah memperlihatkan pengaruh dari pemanenan Predator terhadap dinamik dari sistem. Sistem memiliki paling banyak dua titik tetap yang  memberikan  dinamik  berbeda  disekitar  titik  tetap  ketika  pemanenannya  divariasikan,  dan memperlihatkan  terjadinya bifurkasi  Hopf pada  sistem.  Mortoja et  al .  dalam [5] mengemukakan bahwa Monod dan Haldane melakukan suatu penelitian terkait fungsi respon yang dikenal sebagai fungsi respon Holling tipe IV. Penelitian ini didasari bahwa interaksi predator-prey tidak selamanya monoton,  yaitu  saat  dimana  populasi  mangsa  meningkat  maka  populasi  pemangsa  akan  menurun, dikarenakan  adanya  sifat  bertahan  atau defending dari  mangsa.  Berikut  persamaan  fungsi  respon Monod-Haldane [6] :
𝑝(𝑥) = 𝑎𝑥 𝑛 + 𝛾𝑥 + 𝑥 2
Menurut Ruan dan Xiao, dalam [7] , Sokol dan Howell [8] juga meneliti tentang predator-prey yang  bersifat  tidak  monoton.  Dalam  penelitiannya,  Sokol  dan  Howell  menyatakan  bahwa  model fungsi respon Monod-Haldane yang mereka teliti secara signifikan lebih baik dan lebih sederhana. Berikut persamaannya:
𝑝(𝑥) = 𝑎𝑥 𝑛 + 𝑥 2
Model predator-prey yang melibatkan perlakuan khusus terhadap predator juga telah dibahas pada [9][10]. Dalam artikel ini, kami akan memodifikasi model predator-prey yang ada dalam [1] dengan menggunakan model pertumbuhan logistik [11] , fungsi respon Monod-Haldane dan perilaku anti  pemangsa.  Dari  hasil  modifikasi  model,  kami  melakukan  proses  penentuan  titik  tetap, menganalisis kestabilan dari setiap titik tetap dan melakukan simulasi numerik pada salah satu titik tetap.  Simulasi  tersebut  bertujuan  untuk  memperlihatkan  bahwa  model  memiliki  kestabilan  yang bersifat simpul stabil dan spiral stabil.

## 2. Model Matematika
Model yang akan dibahas dalam artikel ini adalah model predator-prey dengan fungsi respon Monod-Haldane dengan perilaku anti pemangsa. Banyak populasi pada waktu (𝑡) terbagi menjadi 2 populasi yaitu populasi mangsa (𝑥) dan populasi pemangsa (𝑦) . Asumsi-asumsi yang digunakan pada model predator-prey sebagai berikut:
1) Hanya terdapat masing-masing 1 jenis mangsa dan 1 jenis pemangsa dalam populasi.
2) Populasi mangsa diasumsikan bertumbuh secara logistik dan berkurang dengan adanya interaksi antar mangsa dan pemangsa yang tidak monoton.
3) Populasi pemangsa diasumsikan bertambah dengan adanya interaksi antar mangsa dan pemangsa yang tidak monoton, berkurang karena laju kematian alami dan berkurang karena adanya perilaku anti pemangsa dari populasi mangsa.
Berdasarkan asumsi diatas, diperoleh model persamaan berikut:
𝑥̇(𝑡) = 𝑟𝑥 (1 − 𝑥 𝐾 ) − 𝑎𝑥𝑦 𝑥 + 𝑛𝑦 2 𝑦̇(𝑡) = 𝑏𝑥𝑦 𝑥 + 𝑛𝑦 2 − 𝑚𝑦 − 𝜂𝑥𝑦
Keterangan parameter disajikan pada Tabel 1.
(1)
Tabel 1. Parameter yang digunakan
Simbol Keterangan Syarat 𝑥 Banyaknya populasi mangsa 𝑥 ≥ 0 𝑦 Banyaknya populasi pemangsa 𝑦 ≥ 0 𝑟 Laju pertumbuhan intrinsik mangsa 𝑟 ≥ 0 𝐾 Daya tampung lingkungan 𝐾 ≥ 0 𝑎 Koefisien  interaksi  antara  mangsa  dan  pemangsa  yang berpengaruh terhadap laju pertumbuhan mangsa 𝑎 ≥ 0 𝑏 Koefisien  interaksi  antara  mangsa  dan  pemangsa  yang berpengaruh terhadap laju pertumbuhan pemangsa 𝑏 ≥ 0 𝑚 Laju kematian alami pemangsa 𝑚 ≥ 0 𝑛 Tingkat kejenuhan pemangsaan 𝑛 ≥ 0 𝜂 Perilaku anti pemangsa 𝜂 ≥ 0
Dengan menggunakan bentuk simplifikasi, dimana 𝛽 = 1 𝑥+𝑛𝑦 2 , maka persamaan (1) menjadi:
𝑥̇(𝑡) = 𝑟𝑥 (1 − 𝑥 𝐾 ) − 𝑎𝑥𝑦 𝑥 + 𝑛𝑦 2 𝑦̇(𝑡) = 𝑏𝑥𝑦 𝑥 + 𝑛𝑦 2 − 𝑚𝑦 − 𝜂𝑥𝑦
Sesuai  yang  ada  dalam [12] ,  persamaan (2) ditransformasikan  ke  bentuk  yang  lebih  sederhana dengan cara penondimensionalan model [13] . Skala parameter yang digunakan yaitu:
𝑥 → 𝐾𝑥        𝑦 → 𝐾𝑦        𝑡 → 𝑡 𝐾𝑎
Sehingga persamaan (2) menjadi:
𝑥̇(𝑡) = 𝛼𝑥(1 − 𝑥) − 𝛽𝑥𝑦
𝑦̇(𝑡) = 𝛾𝑥𝑦 − 𝛿𝑦 − 𝜙𝑥𝑦
Dengan
𝑎 = 𝑟 𝐾𝑎 𝛾 = 𝛽𝑏 𝑎 𝛿 = 𝑚 𝐾𝑎 𝜙 = 𝜂 𝑎
Dalam persamaan (3) , parameter 𝛼 merepresentasikan tingkat pertumbuhan mangsa, parameter 𝛾 merepresentasikan tingkat interaksi antara mangsa dan pemangsa, parameter 𝛿 merepresentasikan tingkat kematian pemangsa, serta parameter 𝜙 merepresentasikan tingkat perilaku anti pemangsa.

## 3. Hasil dan Pembahasan

## 3.1. Penentuan Titik Tetap
Titik tetap persamaan (3) diperoleh dengan menyelesaikan persamaan 𝑥̇(𝑡) = 0 dan 𝑦̇(𝑡) = 0 , sehingga diperoleh:
𝛼𝑥(1 − 𝑥) − 𝛽𝑥𝑦 = 0
(2)
(3)
(3)

## 𝛾𝑥𝑦 − 𝛿𝑦 − 𝜙𝑥𝑦 = 0
Dengan  menyelesaikan  sistem  persamaan (4) ,  diperoleh  3  titik  tetap,  yaitu 𝑇 1 (0,0) , 𝑇 2 (1,0) dan 𝑇 3 (𝑥 ∗ , 𝑦 ∗ ) , dengan:
𝑥 ∗ = 𝛿 𝛾 − 𝜙 𝑦 ∗ = 𝛼(𝛾 − 𝜙 − 𝛿) 𝛽(𝛾 − 𝜙)
Titik tetap 𝑇 3 (𝑥 ∗ , 𝑦 ∗ ) akan berada dikuadran pertama jika memenuhi:
𝑥 ∗ > 0   atau 𝛿 𝛾 − 𝜙 > 0 𝑦 ∗ > 0   atau 𝛼(𝛾 − 𝜙 − 𝛿) 𝛽(𝛾 − 𝜙) > 0
Agar persamaan (5) dan (6) terpenuhi, maka:
𝛿 > 0   dan  𝛾 > 𝛿 + 𝜙
Pada saat 𝛾 = 𝛿 + 𝜙 , titik tetap 𝑇 3 akan menyatu dengan titik tetap 𝑇 2 . Pada saat 𝛿 = 0 , titik tetap 𝑇 3 akan menyatu dengan titik tetap 𝑇 1 .

## 3.2.  Analisis Kestabilan Titik Tetap
Dengan melakukan pelinearan [14] pada persamaan (4) , diperoleh matriks jacobi sebagai berikut:
𝐽 = ( 𝛼 − 2𝛼𝑥 − 𝛽𝑦 −𝛽𝑥 𝛾𝑦 − 𝜙𝑦 𝛾𝑥 − 𝛿 − 𝜙𝑥 )
Kita dapat menentukan kestabilan dari setiap titik tetap dengan menggunakan persamaan karakteristik [15] dan berdasarkan beberapa lemma berikut [16] :
Lemma 1. Persamaan (4) disekitar titik tetap 𝑇 1 selalu bersifat sadel.
Bukti:
Jika titik tetap 𝑇 1 (0,0) disubstitusikan ke dalam matriks jacobi (7) maka diperoleh:
𝐽 (0,0) = ( 𝛼 0 0 −𝛿 )
Dengan menyelesaikan persamaan karakteristik det(𝐽 (0,0) − 𝜆𝐼) = 0 , diperoleh nilai eigen untuk matriks 𝐽 (0,0) , yaitu:
𝜆 1 = 𝛼  dan  𝜆 2 = −𝛿
Parameter 𝛼 dan 𝛿 diasumsikan bernilai positif, maka 𝜆 1 > 0 dan 𝜆 2 < 0 . Karena nilai-nilai eigen yang diperoleh berbeda tanda, maka titik tetap 𝑇 1 (0,0) memiliki kestabilan yang bersifat sadel. ∎
Lemma 2 . Persamaan (4) disekitar titik tetap 𝑇 2 akan stabil jika 𝛾 < 𝛿 + 𝜙 .
Bukti :
Jika titik tetap 𝑇 2 (1,0) disubstitusikan ke dalam matriks jacobi (7) maka diperoleh:
𝐽 (1,0) = ( −𝛼 −𝛽 0 𝛾 − 𝛿 − 𝜙 )
(5) (6)
(7) (4)
Dengan menyelesaikan persamaan karakteristik det(𝐽 (1,0) − 𝜆𝐼) = 0 , diperoleh nilai eigen untuk matriks 𝐽 (1,0) , yaitu:
𝜆 1 = −𝛼  dan  𝜆 2 = 𝛾 − 𝛿 − 𝜙
Parameter 𝛼 diasumsikan bernilai positif sehingga 𝜆 1 < 0 . Apabila 𝛾 < 𝛿 + 𝜙 , maka diperoleh nilai 𝜆 2 < 0 . Sehingga titik tetap 𝑇 2 bersifat stabil. ∎
Lemma 3. Persamaan (4) disekitar titik tetap 𝑇 3 akan stabil jika 𝐴 < 0 dan 𝐵 > 0 .
Bukti :
Jika titik tetap 𝑇 3 (𝑥 ∗ , 𝑦 ∗ ) disubstitusikan ke dalam matriks jacobi (7) maka diperoleh:
𝐽 (𝑥 ∗ ,𝑦 ∗ ) = ( 𝑀 11 𝑀 12 𝑀 21 𝑀 22 )
Dimana:
𝑀 11 = − 𝛼𝛿 𝛾 − 𝜙 ; 𝑀 12 = − 𝛽𝛿 𝛾 − 𝜙 ; 𝑀 21 = 𝛼(𝛾 − 𝜙 − 𝛿) 𝛽 ; 𝑀 22 = 0
Nilai eigen diperoleh dengan menyelesaikan persamaan karakteristik det(𝐽 (𝑥 ∗ ,𝑦 ∗ ) − 𝜆𝐼) = 0 yaitu:
| 𝑀 11 − 𝜆 𝑀 12 𝑀 21 𝑀 22 − 𝜆 | = 0 𝜆 2 − (𝑀 11 + 𝑀 22 )𝜆 + 𝑀 11 𝑀 22 − 𝑀 12 𝑀 21 = 0 𝜆 2 − 𝐴𝜆 + 𝐵 = 0
Sehingga
𝜆 1,2 = 𝐴 ± √𝐴 2 − 4𝐵 2
Dengan
𝐴 = − 𝛼𝛿 𝛾 − 𝜙 ;    𝐵 = 𝛿𝛼(𝛾 − 𝜙 − 𝛿) 𝛾 − 𝜙
Karena parameter 𝛼 , 𝛿 diasumsikan bernilai positif maka harus memenuhi 𝐴 < 0 . Jika 𝐵 > 0 , diperoleh 𝛿 < 𝛾 − 𝜙 . Hal ini sesuai dengan syarat (21) , maka titik tetap 𝑇 3 (𝑥 ∗ , 𝑦 ∗ ) bersifat stabil. ∎

## 3.3.  Simulasi Numerik
Dinamika populasi mangsa-pemangsa dapat ditunjukkan melalui kurva bidang fase dan bidang solusi yang menggambarkan populasi mangsa dan pemangsa pada kurun waktu tertentu. Dengan ini, masing-masing  variabel  dan  parameter  membutuhkan  suatu  nilai  awal  untuk  proses  komputasi. Proses komputasi pada simulasi ini dilakukan dengan bantuan software Wolfram Mathematica 11.3 dan software Maple 18 .
Diasumsikan nilai parameter yang bernilai tetap untuk simulasi ini adalah 𝛼 = 0.1 , 𝛿 = 0.01 dan 𝜙 = 0.2 . Parameter 𝛽 dan 𝛾 merupakan parameter kontrol pada sistem yang nilainya tidak selalu tetap. Parameter yang digunakan dalam simulasi ini dapat dilihat pada Tabel 2.
Tabel 2. Pemilihan Nilai Parameter Parameter Kasus 1 Kasus 2 𝛼 0.1 0.1
𝛿 0.01 0.01 𝜙 0.2 0.2 𝛽 0.22 0.23 𝛾 0.22 0.23 𝑇 3 Simpul Stabil Spiral Stabil

## 3.3.1. Dinamika Populasi Predator-Prey Kasus 1
Pada kasus ini, parameter 56amper56 yang digunakan yaitu 𝛽 = 𝛾 = 0.22 , serta nilai awal 𝑥(0) = 0.3 dan 𝑦(0) = 0.5 . Titik tetap 𝑇 3 yang diperoleh pada kasus ini adalah 𝑇 3 (0.5,0.227273) . Nilai eigen yang  diperoleh  adalah 𝜆 1 = −0.0361803 dan 𝜆 2 = −0.0138197 ,  sehingga  kestabilan  titik  tetap 𝑇 3 bersifat simpul stabil.
Pada Gambar 1, diberikan ilustrasi bidang fase disekitar titik tetap dimana kedua populasi stabil menuju titik tetap 𝑇 3 . Diperlihatkan bahwa jenis kestabilan titik tetapnya adalah stabil.
Pada Gambar 2, diperlihatkan bahwa di awal waktu kedua populasi mengalami penurunan jumlah, dimana populasi  mangsa  56amper  mengalami  kepunahan  dikarenakan  sifat defending yang  dimiliki.  Setelah  itu,
populasi pemangsa mengalami penurunan yang mengakibatkan populasi mangsa mengalami pertumbuhan. Pada akhirnya kedua populasi tersebut stabil menuju ke satu nilai.

## 3.3.2. Dinamika Populasi Predator-Prey Kasus 2
Pada kasus ini, parameter kontrol yang digunakan yaitu 𝛽 = 𝛾 = 0.23 , serta nilai awal 𝑥(0) = 0.2 dan 𝑦(0) = 0.6 . Titik tetap 𝑇 3 yang diperoleh pada kasus ini adalah 𝑇 3 (0.333333,0.289855) . Nilai eigen yang  diperoleh  adalah 𝜆 1 = −0.0166667 + 0.0197203𝒊 dan 𝜆 2 = −0.0166667 + 0.0197203𝒊 , sehingga kestabilan titik tetap 𝑇 3 bersifat spiral stabil.
Pada Gambar 3, diberikan ilustrasi bidang fase disekitar titik tetap dimana kedua populasi stabil menuju titik tetap 𝑇 3 . Diperlihatkan bahwa jenis kestabilan titik tetapnya adalah stabil.
Pada Gambar 4, diperlihatkan bahwa di awal waktu kedua populasi mengalami penurunan jumlah, dimana populasi mangsa hampir mengalami kepunahan pada interval waktu 0 − 50 hari dikarenakan sifat defending yang dimiliki sedangkan populasi pemangsa mengalami penurunan jumlah pada interval waktu 0 − 150 hari. Kemudian pada interval waktu 50 − 200 hari populasi mengalami peningkatan jumlah populasi yang  signifikan  sedangkan  pada  interval  waktu 150 − 300 hari  populasi  pemangsa  juga  mengalami peningkatan dikarenakan jumlah populasi mangsa meningkat sehingga menimbulkan interaksi antar mangsa
dan pemangsa. Oleh sebab itu, pada interval waktu 200 − 300 hari populasi mangsa mengalami penurunan dikarenakan interaksi tersebut. Pada akhirnya kedua populasi tersebut stabil menuju ke satu nilai.

## 4. Kesimpulan
Dari hasil analisis model, diperoleh tiga titik tetap, yaitu 𝑇 1 , 𝑇 2 dan 𝑇 3 . Kestabilan titik tetap 𝑇 1 selalu  bersifat  sadel.  Kestabilan  titik  tetap 𝑇 2 bergantung  pada  nilai  parameter  yang  dipilih  yang menyebabkan  kestabilannya  bersifat  sadel  jika 𝛾 > 𝛿 + 𝜙 dan  bersifat  stabil  jika 𝛾 < 𝛿 + 𝜙 . Kestabilan  titik  tetap 𝑇 3 juga  akan  stabil  jika  memenuhi  syarat 𝐴 < 0 dan 𝛿 < 𝛾 − 𝜙 .  Dinamika populasi mangsa-pemangsa model ini terbagi atas 2 kasus, dimana parameter 𝛽 dan 𝛾 merupakan parameter kontrol pada sistem yang nilainya tidak selalu tetap. Pada kasus 1, kestabilan titik tetap 𝑇 3 bersifat simpul stabil dimana nilai parameter kontrolnya yaitu 𝛽 = 𝛾 = 0.22 . Sedangkan pada kasus 2, kestabilan titik tetap 𝑇 3 bersifat spiral stabil dimana nilai parameter kontrolnya yaitu 𝛽 = 𝛾 = 0.23 . Pada kedua kasus ini, hampir terjadi kepunahan pada populasi mangsa yang disebabkan oleh perilaku anti pemangsa dan sifat defending dari mangsa.

## Referensi
[1] S.  Saadah,  A.  Abadi,  and  D.  Savitri, “ Model  Interaksi  Mangsa  Pemangsa  Dengan  Fungsi Respon Rasio Dependent Holling Tipe II dan Perilaku Anti Pemangsa, ” MATHunesa Jurnal Ilmiah Matematika, vol.7, no.2, pp. 6-9, 2019.
[2] P.B. Turchin, Complex Population Dynamics: A Theorical/Empirical Synthesis. New Jersey: Princeton University Press, 2003.
[3] W. Widowati dan S. Sutimin, Buku Ajar Pemodelan Matematika, Semarang: UNDIP  Press, 2007.
[4] H.S.  Panigoro, “ Analisis  Dinamik  Sistem  Predator-Prey  Model  Leslie-Gower  dengan Pemanenan Secara Konstan Terhadap Predator, ” J.Euler, vol.2, no.1, 1-12, 2014.
[5] S.G. Mortoja, P. Panja, and S.K. Mondal, “ Dynamics of a Predator-Prey Model with Stage- Structure  on  Both  Species  and  Anti-Predator  Behavior, ” Informatics  in  Medicine Unlocked , vol.10 , pp. 50-57, 2018. doi:https://doi.org/10.1016/j.imu.2017.12.004.
[6] M. Agarwal, and R. Pathak, “ Harvesting and Hopf Bifurcation in a Predator-Prey Model with Holling  Type  IV  Functional  Response, ” International  Journal  of  Mathematics  and  Soft Computing, vol.2, no.1, pp. 83-92, 2012.
[7] S.  Ruan,  and  D.  Xiao, “ Global  Analysis  in  a  Predator-Prey  System  with  Nonmonotonic Functional Response, ” SIAM J Appl Math , vol. 61, no.4, pp. 1445-1472, 2001.
[8] W. Sokol, and J.A. Howell, “ Kinetics of Phenol Oxidation by Washed Cells, ” Biotechnol . Bio- eng. , vol. 23, pp. 2039-2049, 1980.
[9] S. H. A rsyad, R. Resmawan, and N. Achmad, “Analisis Model Predator -Prey Leslie-Gower dengan Pemberian Racun Pada Predator,” J. Ris. dan Apl. Mat. , vol. 4, no. 1, pp. 1 – 16, 2020.
[10] N. Hasan, R. Resmawan, and E. Rahmi, “Analisis Kestabilan Model Eko -Epidemiologi dengan Pemanenan Konstan pada Predator,” J. Mat. Stat. dan Komputasi , vol. 16, no. 2, pp. 121 – 142, Dec. 2020.
[11] J.D.  Logan, First  Course  in  Differential  Equations  (Second  Edition), New  York:  Springer- Verlag, 2010.
[12] M.B.  Gaib, Bifurkasi  Hopf  Pada  Model Mangsa-Pemangsa  Holling-Tanner  Tipe  II .  Institut Pertanian Bogor, 2013.
[13] S.H. Strogatz, Nonlinear Dynamics and Chaos with Application to Physics, Biology, Chemistry, and Engineering , Massachusets (US): Addison-Wesley Publishing Company, 1994.
[14] J.K. Hale, and H. Kocak, Dynamic and Bifurcation, New York: Springer-Verlag, 1991.
[15] H. Anton, and C. Rorres, Aljabar Linear Elementer, Jakarta: Erlangga, 1987.
[16] M.  Diska, Analisis  Kestabilan  Model  Mangsa-Pemangsa  Dengan  Struktur  Umur  Pemangsa dan Fungsi Respon Monod-Haldane, Institut Pertanian Bogor, 2018.