## RESEARCH ARTICLE • OPEN ACCESS
Efek Perlindungan Mangsa dan Daya Dukung Variabel pada Sistem Mangsa-Pemangsa dengan Fungsi Respon Beddington-DeAngelis

## Aisyiah Kholifatul Hanifah dan Abadi Abadi

## Volume 6, Issue 1, Pages 102–110, February 2024
Submit 15 Januari 2024, Direvisi 7 Februari 2024, Disetujui 10 Februari 2024 To Cite this Article : A. K. Hanifah dan A. Abadi,“Efek Perlindungan Mangsa dan Daya Du- kung Variabel pada Sistem Mangsa-Pemangsa dengan Fungsi Respon Beddington-DeAngelis ”, Jambura J. Math , vol. 6, no. 1, pp. 102–110, 2024, https://doi.org/10.37905/jjom.v6i1.24288
© 2024 by author(s)

## JOURNAL INFO • JAMBURA JOURNAL OF MATHEMATICS
u Homepage : http://ejurnal.ung.ac.id/index.php/jjom/index Ǽ Journal Abbreviation : Jambura J. Math. Х Frequency : Biannual (February and August) Ŷ Publication Language : English (preferable), Indonesia  DOI : https://doi.org/10.37905/jjom ʢ Online ISSN : 2656-1344 ǽ Editor-in-Chief : Hasan S. Panigoro Ÿ Publisher : Department of Mathematics, Universitas Negeri Gorontalo ƒ Country : Indonesia ̰ OAI Address : http://ejurnal.ung.ac.id/index.php/jjom/oai  Google Scholar ID : iWLjgaUAAAAJ Ͳ Email : info.jjom@ung.ac.id

## JAMBURA JOURNAL • FIND OUR OTHER JOURNALS
Jambura Journal of
Biomathematics Jambura Journal of Mathematics Education Jambura Journal of Probability and Statistics EULER : Jurnal Ilmiah Matematika, Sains, dan Teknologi

## Efek Perlindungan Mangsa dan Daya Dukung Variabel pada Sistem Mangsa-Pemangsa dengan Fungsi Respon

## Beddington-DeAngelis
Aisyiah Kholifatul Hanifah 1 dan Abadi Abadi 1 , ∗ 
1 Jurusan Matematika, Universitas Negeri Surabaya, Indonesia
ABSTRAK. Beberapa spesies di dunia telah mengalami kepunahan. Untuk menyelamatkan spesies dari kepunahan, formulasi dari sistem perlu dilakukan, salah satunya dengan perlindungan mangsa. Perubahan proses fisik dan biolo- gis juga memiliki peran pada lingkungan, yaitu mengakibatkan daya dukung yang berubah secara dinamis. Tingkat konsumsi mangsa oleh pemangsa rata-rata, yaitu fungsi respon juga penting dalam formulasi sistem. Pada penelitian ini dikonstruksikan sistem mangsa-pemangsa dengan fungsi respon Beddington-DeAngelis yang mempertimbangkan perlindungan mangsa dan daya dukung variabel serta proses predasi yang memperhatikan populasi pemangsa. Da- lam penelitian ini ditentukan titik kesetimbangan kepunahan mangsa-pemangsa, kepunahan pemangsa, dan hidup berdampingan lalu dilanjutkan menganalisis kestabilannya. Dengan nilai perlindungan mangsa yang besar, titik kesetimbangan kepunahan pemangsa bersifat stabil asimtotik dan titik kesetimbangan hidup berdampingan bersifat stabil asimtotik, jika dan hanya jika pada rentang nilai parameter perlindungan mangsa tertentu. Selain itu, dila- kukan simulasi dan didapatkan kesimpulan bahwa pada nilai perlindungan mangsa tertentu, terjadi perubahan pada kestabilan sistem yang bergantung pada batas-batas daya dukungnya.
ABSTRACT. Several species in the world have experienced extinction. To save species from extinction, a system needs to be formulated, one of which is protecting prey. Changes in physical and biological processes also have a role in the environment, namely resulting in dynamically changing carrying capacity. The rate of prey consumption by the average predator, i.e. the response function is also important in the system formulation. In this study, a prey- predator system was constructed with a Beddington-DeAngelis response function that considers prey protection and variable carrying capacity as well as a predation process that takes into account the predator population. In this research, the equilibrium point for prey-predator extinction, predator extinction, and coexistence was determined and then continued to analyze its stability. With a large prey protection value, the predator extinction equilibrium point is asymptotically stable and the coexistence equilibrium point is asymptotically stable, if and only if within a certain range of prey protection parameter values. In addition, simulations were carried out and it was concluded that at certain prey protection values, changes occurred in the stability of the system which depended on the limits of its carrying capacity.
This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution-NonComercial 4.0 International License. Editorial of JJBM: Department of Mathematics, Uni- versitas Negeri Gorontalo, Jln. Prof. Dr. Ing. B. J. Habibie, Bone Bolango 96554, Indonesia.

## 1. Pendahuluan
Model matematika sering kali membahas permasalahan fe- nomena alam yang terjadi pada makhluk hidup, serta dapat di- gunakan untuk memprediksi perilaku sistem untuk jangka waktu tertentu yang direpresentasikan dalam suatu sistem persamaan. Salah satu model interaksi makhluk hidup dalam satu ekosistem adalah model mangsa-pemangsa.
Beberapa spesies di dunia ini telah mengalami kepunahan dan banyak spesies lainnya yang hampir punah. Hal ini disebabk- an oleh berbagai alasan misalnya pencemaran lingkungan, eks- ploitasi, pemanenan yang tidak teratur, pencemaran lingkungan, tingkat pemangsaan yang terlalu tinggi. Untuk menyelamatkan
∗ Penulis Korespondensi.
Check for updates
Research Article

## ARTICLE HISTORY
Submit 15 Januari 2024
Direvisi 7 Februari 2024 Disetujui 10 Februari 2024

## KATA KUNCI
Sistem Mangsa-Pemangsa
Beddington-DeAngelis Analisis Kestabilan
Perlindungan Mangsa
Daya dukung Variabel

## KEYWORDS
Prey-Predator System Beddington-DeAngelis
Stability Analysis
Prey Protection
Variable Carrying Capacity
spesies dari kepunahan, formulasi dari sistem perlu dilakukan [1– 3] .
Perilaku pada mangsa salah satunya yaitu melindungi diri dari pemangsa atau dapat disebut sebagai perlindungan mang- sa. Beberapa studi terdahulu telah menunjukkan efek dari men- stabilkan dinamika sistem dengan suatu perlindungan mangsa. Contohnya adalah penelitian yang dilakukan oleh Connell dkk pada tahun 1970 [4] , ditemukan fakta bahwa perlindungan spa- sial yang ditunjukkan oleh spesies teritip ( Balanus glandula, Ba- lanus cariosus, Chthamalus fissus ) di daerah pasang surut yang le- bih tinggi dapat menstabilkan interaksinya dengan siput preda- tor ( Thais emarginata, Thais canaliculata, Thais lamellosa ). Seca- ra umum, perlindungan mangsa memiliki peran penting dalam kestabilan interaksi mangsa-pemangsa dan dapat mencegah ter-
Email : abadi@unesa.ac.id (A. Abadi) Homepage : http://ejurnal.ung.ac.id/index.php/jjom/index / E-ISSN : 2656-1344 © 2024 by the Author(s).
jadinya kepunahan spesies mangsa, karena mangsa melindungi diri dari pemangsa. Beberapa penulis telah mempelajari dam- pak tempat perlindungan mangsa terhadap stabilitas dinamika mangsa-pemangsa [2, 3, 5–11] .
Selain perlindungan mangsa yang memiliki peran penting dalam kestabilan sistem mangsa-pemangsa, perubahan proses fisik dan biologis juga memiliki peran pada lingkungan, yaitu mengakibatkan perubahan daya dukung dan pada akhirnya mem- berikan pengaruh pada sistem mangsa-pemangsa. Di sebagian besar literatur atau penelitian terdahulu misalnya [1, 2, 12] , daya dukung dinyatakan dalam bentuk konstan. Tetapi, dalam keada- an lingkungan yang nyata, daya dukung berubah secara dinamis, jadi masuk akal jika mempertimbangkan daya dukung sebagai va- riabel model (yaitu fungsi dalam waktu). Selain itu, pada kenyata- annya sumber daya bersifat terbatas sehingga peningkatan atau penurunan daya dukung memiliki suatu batas [5, 13–15] .
Fungsi respon berkaitan erat dengan peningkatan populasi pemangsa atau pengurangan populasi mangsa saat saling berin- teraksi [16] . Terdapat keadaan di mana peningkatan kepadatan pemangsa yang mengisyaratkan penurunan laju makan pemang- sa karena adanya gangguan timbal balik di antara pemangsa- pemangsa. Pada penelitian yang dilakukan oleh Hassel pada ta- hun 1971 [17] menunjukkan bahwa efisiensi pencarian parasit se- rangga, Nemeritis canescens , menurun pada kepadatan parasit atau spesies serangga hama terutama ulat yang tinggi. Hal ini meru- pakan faktor kunci untuk dapat memodifikasi fungsi respon Ho- lling tipe-II ke dalam bentuk fungsi respon Beddington-DeAngelis . Fungsi respon Beddington-DeAngelis serupa dengan respon fung- sional Holling tipe-II tetapi mengandung istilah tambahan yang menggambarkan gangguan timbal balik oleh pemangsa [18, 19] .
Pada penelitian terdahulu yang telah dilakukan oleh Al-Salti dkk pada tahun 2021 [5] mengenai efek penggabungan perlin- dungan mangsa m dan daya dukung variabel yang dibatasi antara K 1 dan K 2 dalam sistem dinamika interaksi mangsa-pemangsa de- ngan fungsi respon Holling tipe-II, ditemukan bahwa untuk nilai batas awal daya dukung yang kecil, populasi mangsa-pemangsa untuk kasus daya dukung yang menurun mencapai kesetimbang- an kepunahan pemangsa untuk nilai perlindungan mangsa yang tinggi. Dengan meningkatkan batas atas daya dukung, rentang nilai perlindungan mangsa dan batas awal daya dukung membuat sistem mencapai kesetimbangan hidup berdampingan yang sta- bil. Dengan simulasi numerik ditunjukkan bahwa kasus daya du- kung variabel dan daya dukung konstan, memiliki perilaku jangka panjang yang sama namun, dinamika jangka pendeknya berbeda. Selain itu, ditunjukkan juga bahwa ketika daya dukung berubah secara perlahan, sistem membutuhkan waktu yang lebih lama un- tuk mencapai titik kesetimbangan.
Berdasarkan uraian di atas, pada penelitian ini dikaji lebih lanjut mengenai model mangsa-pemangsa yang mengacu pada penelitian Al-Salti dkk [5] dengan menambahkan asumsi bahwa terkait laju pemangsaan tidak hanya bergantung pada mangsa melainkan juga bergantung pada pemangsa. Dalam hal ini dia- sumsikan terjadi peningkatan pada kepadatan populasi pemang- sa yang mengakibatkan penurunan laju pemangsaan sebagaima- na yang dibahas pada penelitian Tripathi dkk [1] . Modifikasi di- lakukan untuk mengetahui bagaimana efek perlindungan mang- sa terhadap fenomena kepunahan populasi mangsa dan daya du- kung variabel pada sistem mangsa-pemangsa dengan fungsi res-
pon Beddington-DeAngelis . Sistem mangsa-pemangsa disajikan ke dalam 3 populasi yaitu N ( Prey /Mangsa), P ( Predator /Pemangsa), dan K ( Carrying Capacity /Daya Dukung).
Penelitian ini dilakukan dengan mengkaji studi literatur, mengkonstruksi model mangsa-pemangsa, menentukan titik ke- setimbangan, melakukan analisis kestabilan dan analisis bifurka- si, serta membuat simulasi numerik untuk mengonfirmasi hasil analisis yang telah dilakukan.

## 2. Model
Laju populasi mangsa N terhadap waktu dipengaruhi oleh populasi mangsa yang tumbuh secara logistik dengan tingkat per- tumbuhan r dengan kata lain populasi mangsa tumbuh secara eks- ponensial sebesar rN dan akan berkurang karena adanya faktor daya dukung N 2 K . Selain itu populasi mangsa akan berkurang dengan tingkat penurunan mangsa a akibat dari proses preda- si yang mengalami hambatan karena adanya faktor penanganan pemangsa saat melakukan pemangsaan γ 1 dan adanya gangguan antar populasi pemangsa γ 2 . Selain itu, jumlah populasi mang- sa yang dimangsa tersisa (1 − m ) n karena adanya faktor perlin- dungan mangsa. Dengan demikian laju populasi mangsa atau N terhadap waktu dapat direpresentasikan pada persamaan (1) :
dN dt = rN ( 1 − N K ) − a (1 − m ) N P 1 + γ 1 (1 − m ) N + γ 2 P . (1)
Laju populasi pemangsa P terhadap waktu akan bertambah dengan tingkat kenaikan pemangsa b dipengaruhi oleh proses predasi yang mengikuti fungsi respon Beddington-DeAngelis di ma- na selain memperhatikan populasi mangsa juga memperhatikan gangguan di antara pemangsa. Populasi pemangsa memakan po- pulasi mangsa untuk mendapat energi dimana yang selain digu- nakan untuk hidup, energi tersebut digunakan juga untuk berep- roduksi menghasilkan populasi pemangsa baru. Selain itu popu- lasi pemangsa akan berkurang akibat dari proses kematian alami pemangsa dengan tingkat kematian alami pemangsa c , sehingga laju populasi pemangsa atau P terhadap waktu dapat direpresen- tasikan pada persamaan (2) :
dP dt = b (1 − m ) N P 1 + γ 1 (1 − m ) N + γ 2 P − cP. (2)
Daya dukung diasumsikan berubah secara dinamis mengi- kuti model pertumbuhan logistik tergantung pada nilai awal dan nilai pembatasnya K 1 , K 2 . Daya dukung berubah dengan ting- kat pertumbuhan α dengan kata lain daya dukung tumbuh se- cara eksponensial dengan batas bawah daya dukung K 1 yaitu α ( K − K 1 ) dan akan berkurang karena adanya batas atas dan batas bawah daya dukung α ( K − K 1 ) 2 δ k ,di mana δ K = K 2 – K 1 . De- ngan demikian laju daya dukung atau K terhadap waktu dapat dapat direpresentasikan pada persamaan (3) :
dK dt = α ( K − K 1 ) ( 1 − K − K 1 δ k ) . (3)
Dari persamaan (1) - (3) , diperoleh sistem mangsa-pemangsa dengan fungsi respon Beddington-DeAngelis yang mempertim- bangkan perlindungan mangsa dan daya dukung variabel pada persamaan (4) :
dN dt = rN ( 1 − N K ) − a (1 − m ) N P 1 + γ 1 (1 − m ) N + γ 2 P dP dt = b (1 − m ) N P 1 + γ 1 (1 − m ) N + γ 2 P − cP dK dt = α ( K − K 1 ) ( 1 − K − K 1 δ K ) . (4)
dengan keterangan variabel atau parameter yang digunakan dije- laskan pada Tabel 1 .
Tabel 1. Keterangan variabel dan parameter Variabel Keterangan N Populasi mangsa P Populasi pemangsa K Daya dukung r Laju pertumbuhan mangsa a Tingkat penurunan mangsa b
Tingkat kenaikan pemangsa
c Tingkat kematian alami pemangsa
m Konstanta perlindungan mangsa, m ∈ [0, 1) γ 1 Parameter penanganan γ 2 Besarnya gangguan di antara pemangsa α Tingkat penurunan atau kenaikan daya dukung
K 1 , K 2 Batas daya dukung, δ K = K 2 − K 1 Sumber: [1, 5]
Untuk mempermudah penentuan titik kesetimbangan, jum- lah parameter dalam sistem atau model (4) dapat direduksi de- ngan penskalaan mengikuti persamaan (5) :
n = N n c , p = P p c , k = K k c , τ = t τ c . (5)
Sistem (4) dapat ditransormasi menggunakan penskalaan (5) , sehingga diperoleh sistem berskala pada persamaan (6) :
dn dτ = n ( 1 − n K ) − (1 − m ) np 1 + β 1 (1 − m ) n + β 2 p dp dτ = (1 − m ) np 1 + β 1 (1 − m ) n + β 2 p − σ 1 p dk dτ = σ 2 ( k − k 1 ) ( 1 − k − k 1 δ k ) (6)
dimana,
n c = k c = r b , p c = r a , τ c = 1 r
dan
β 1 = rγ 1 b , β 2 = rγ 2 a , σ 1 = c r , σ 2 = α r δ k = k 2 − k 1 , k 1 = rK 1 b , k 2 = rK 2 b .

## 3. Hasil dan Pembahasan
Pada bagian ini, disajikan keberadaan titik keseimbangan sistem (6) , analisis titik kesetimbangannya, serta simulasi nume- riknya.

## 3.1. Titik Kesetimbangan
Titik kesetimbangan diperoleh dengan menyelesaikan sis- tem (6) , yaitu kondisi yang diperoleh ketika dn dτ = dp dτ = dk dτ = 0 . Terdapat 6 kemungkinan titik kesetimbangan yang dapat dike- lompokkan dalam 3 jenis titik kesetimbangan sebagai berikut: 1. Titik kesetimbangan kepunahan mangsa-pemangsa
ξ 1 = (0 , 0 , k 1 ) dan ξ 2 = (0 , 0 , k 2 )

## 2. Titik kesetimbangan kepunahan pemangsa
ξ 3 = ( k 1 , 0 , k 1 ) dan ξ 4 = ( k 2 , 0 , k 2 )

## 3. Titik kesetimbangan hidup berdampingan
ξ 5 = ( n ∗ 1 , p ∗ 1 , k 1 ) dan ξ 6 = ( n ∗ 2
, p ∗ 2 , k 2 )
dimana n ∗ i , i = 1 , 2 adalah solusi dari persamaan kuadrat (7) :
An 2 + Bn + C = 0 . (7)
dengan
A = β 2 , B = β 1 k i mσ 1 − β 1 k i σ 1 − β 2 k i − k i m + k i , i = 1 , 2 C = − k i σ 1 , i = 1 , 2 .
dan
p ∗ i = (1 − m ) (1 − σ 1 β 1 ) n ∗ i − σ 1 σ 1 β 2 , i = 1 , 2 .
Titik kesetimbangan ξ 5 dan ξ 6 ada ketika
(1 − m ) > σ 1 (1 − β 1 σ 1 ) n ∗ , i i = 1 , 2 .

## 3.2. Analisis Kestabilan
Kestabilan titik–titik kesetimbangan ditentukan berdasark- an nilai eigen dari matriks Jacobian pada setiap titik kesetimbang- an sistem (6) . Untuk tujuan ini, kita perhatikan bahwa matriks Jacobian sistem (6) pada titik kesetimbangan, dilambangkan de- ngan J ( ξ ) , yaitu
J ( ξ ) =   J 11 J 12 J 13 J 21 J 22 0 0 0 J 33   (8)
dimana
J 11 = 1 − 2 n k − ( (1 − m ) p + β 2 (1 − m ) p 2 (1 + β 1 (1 − m ) n + β 2 p ) 2 ) J 12 = − ( (1 − m ) n + β 1 (1 − m ) 2 n 2 (1 + β 1 (1 − m ) n + β 2 p ) 2 ) J 13 = n 2 k 2 J 21 = (1 − m ) p + β 2 (1 − m ) p 2 (1 + β 1 (1 − m ) n + β 2 p ) 2 J 22 = (1 − m ) n + β 1 (1 − m ) 2 n 2 (1 + β 1 (1 − m ) n + β 2 p ) 2 − σ 1 J 33 = σ 2 ( 1 − k − k 1 δ k ) − σ 2 ( k − k 1 δ k ) .
Berdasarkan matriks Jacobian pada persamaan (8) , dilakuk- an analisis kestabilan pada masing-masing titik kesetimbangan sebagai berikut:

## 3.2.1. Titik Kesetimbangan Kepunahan Mangsa-Pemangsa
Sifat kestabilan untuk titik kesetimbangan kepunahan mangsa-pemangsa dinyatakan dalam Teorema 1 .
Teorema 1. Titik kesetimbangan kepunahan mangsa-pemangsa ξ 1 = (0 , 0 , k 1 ) dan ξ 2 = (0 , 0 , k 2 ) merupakan titik pelana yang tidak stabil.
Bukti. Dengan mengevaluasi matriks Jacobian (8) di titik kesetim- bangan ξ 1 dan ξ 2 diperoleh:
1. Matriks Jacobian (8) di titik ξ 1 , yaitu
J ( ξ 1 ) =   1 0 0 0 − σ 1 0 0 0
σ 2   . (9)
Dari matriks Jacobian persamaan (9) , diperoleh nilai eigen masing-masing λ 1 = 1 , λ 2 = − σ 1 , λ 3 = σ 2 , sehingga titik kesetimbangan ξ 1 merupakan titik pelana yang tidak stabil karena λ 1 = 1 dan λ 3 = σ 2 positif. 2. Matriks Jacobian (8) di titik ξ 2 , yaitu
J ( ξ 2 ) =   1 0 0 0 − σ 1 0 0 0 − σ 2   . (10)
Dari matriks Jacobian persamaan (10) , diperoleh nilai eigen masing-masing λ 1 = 1 , λ 2 = − σ 1 , λ 3 = − σ 2 , sehingga titik kesetimbangan ξ 2 merupakan titik pelana yang tidak stabil karena λ 1 = 1 positif.

## 3.2.2. Titik Kesetimbangan Kepunahan Pemangsa
Sifat kestabilan untuk titik kesetimbangan kepunahan pe- mangsa dinyatakan dalam Teorema 2 .
Teorema 2. Titik kesetimbangan kepunahan pemangsa ξ 3 = ( k 1 , 0 , k 1 ) selalu tidak stabil dan merupakan titik pelana, se- dangkan titik kesetimbangan ξ 4 = ( k 2 , 0 , k 2 ) stabil asimtotik lokal dan merupakan titik simpul jika m > 1 − σ 1 k 2 (1 − σ 1 β 1 ) .
Bukti. Dengan mengevaluasi matriks Jacobian (8) di titik kesetim- bangan ξ 3 dan ξ 4 diperoleh:
1. Matriks Jacobian (8) di titik ξ 3 , yaitu
J ( ξ 3 ) =    − 1 − (1 − m ) k 1 1+ β 1 (1 − m ) k 1 1 0 (1 − m ) k 1 1+ β 1 (1 − m ) k 1 − σ 1 0
0 0 σ 2    . (11)
Dari matriks Jacobian persamaan (11) , diperoleh nilai eigen
masing-masing:
λ 1 = − 1 ,
λ 2 = ((1 − m ) k 1 + β 1 (1 − m ) 2 k 1 2 ) (1 + β 1 (1 − m ) k 1 ) 2 − σ 1 , λ 3 = σ 2 .
sehingga titik kesetimbangan ξ 3 merupakan titik pelana yang tidak stabil karena λ 3 = σ 2 positif.
2. Matriks Jacobian (8) di titik ξ 4 , yaitu
J ( ξ 4 ) =    − 1 − (1 − m ) k 2 1+ β 1 (1 − m ) k 2 1 0 (1 − m ) k 2 1+ β 1 (1 − m ) k 2 − σ 1 0 0 0 − σ 2    . (12)
Dari matriks Jacobian persamaan (12) , diperoleh nilai eigen masing-masing:
λ 1 = − 1 ,
λ 2 = ((1 − m ) k 1 + β 1 (1 − m ) 2 k 1 2 ) (1 + β 1 (1 − m ) k 1 ) 2 − σ 1 , λ 3 = − σ 2
sehingga titik kesetimbangan ξ 4 merupakan titik simpul yang stabil asimtotik lokal jika λ 2 < 0 atau m > 1 − σ 1 k 2 (1 − σ 1 β 1 ) .

## 3.2.3. Titik Kesetimbangan Hidup Berdampingan
Sifat kestabilan untuk titik kesetimbangan hidup berdam- pingan antara kedua populasi dinyatakan dalam Teorema 3 .
Teorema 3. Titik kesetimbangan hidup berdampingan ξ 5 = ( n ∗ 1 , p ∗ 1 , k 1 ) tidak stabil, sedangkan titik kesetimbangan ξ 6 = ( n ∗ 2 , p ∗ 2 , k 2 ) stabil asimtotik lokal jika trace ( J ∗ ( ξ 6 )) < 0 dan det ( J ∗ ( ξ 6 )) > 0 .
Bukti. Dengan mengevaluasi matriks Jacobian (8) di titik kesetim- bangan ξ 5 dan ξ 6 diperoleh:
1. Matriks Jacobian (8) di titik ξ 5 , yaitu
J ( ξ 5 ) =   J 11 J 12 J 13 J 21 J 22 0 0 0 J 33   , (13)
dimana
J 11 = 1 − 2 n ∗ 1 k 1 − ( (1 − m ) p ∗ 1 + β 2 (1 − m ) p ∗ 1 2 (1 + β 1 (1 − m ) n ∗ 1 + β 2 p ∗ 1 ) 2 )
J 12 = − ( (1 − m ) n ∗ 1 + β 1 (1 − m ) 2 n ∗ 1 2 (1 + β 1 (1 − m ) n ∗ 1 + β 2 p ∗ 1 ) 2 )
J 13 = n ∗ 1 2 k 1 2
J 21 = (1 − m ) p ∗ 1 + β 2 (1 − m ) p ∗ 1 2 (1 + β 1 (1 − m ) n ∗ 1 + β 2 p ∗ 1 ) 2
J 22 = (1 − m ) n ∗ 1 + β 1 (1 − m ) 2 n ∗ 1 2 (1 + β 1 (1 − m ) n ∗ 1 + β 2 p ∗ 1 ) 2 − σ 1 J 33 = σ 2 .
Dari matriks Jacobian persamaan (13) , diperoleh salah satu nilai eigen λ 3 = σ 2 . Karena semua nilai parameter positif, maka λ 3 = σ 2 > 0 . Hal ini sudah cukup untuk menunjukkan bahwa titik kesetimbangan ξ 5 tidak stabil.
2. Matriks Jacobian (8) di titik ξ 6 , yaitu
J ( ξ 6 ) =   J 11 J 12 J 13 J 21 J 22 0 0 0 J 33   , (14)
dimana
J 11 = 1 − 2 n ∗ 2 k 2 − ( (1 − m ) p ∗ 2 + β 2 (1 − m ) p ∗ 2 2 (1 + β 1 (1 − m ) n ∗ 2 + β 2 p ∗ 2 ) 2 ) J 12 = − ( (1 − m ) n ∗ 2 + β 1 (1 − m ) 2 n ∗ 2 2 (1 + β 1 (1 − m ) n ∗ 2 + β 2 p ∗ 2 ) 2 ) J 13 = n ∗ 2 2 k 2 2
J 21 = (1 − m ) p ∗ 2 + β 2 (1 − m ) p ∗ 2 2 (1 + β 1 (1 − m ) n ∗ 2 + β 2 p ∗
2 ) 2
J 22 = (1 − m ) n ∗ 2 + β 1 (1 − m ) 2 n ∗ 2 2 (1 + β 1 (1 − m ) n ∗ 2 + β 2 p ∗ 2 ) 2 − σ 1
J 33 = − σ 2 .
Dari matriks Jacobian persamaan (14) , diperoleh salah satu nilai eigen negatif, yaitu λ 3 = − σ 2 . Dua nilai eigen lain- nya adalah nilai eigen dari matriks tereduksi pada persama- an (15) , yaitu:
J ∗ ( ξ 6 ) = [ J 11 J 12 J 21 J 22 ] , (15)
dimana
J 11 = 1 − 2 n ∗ 2 k 2 − ( (1 − m ) p ∗ 2 + β 2 (1 − m ) p ∗ 2 2 (1 + β 1 (1 − m ) n ∗ 2 + β 2 p ∗ 2 ) 2 ) J 12 = − ( (1 − m ) n ∗ 2 + β 1 (1 − m ) 2 n ∗ 2 2 (1 + β 1 (1 − m ) n ∗ 2 + β 2 p ∗ 2 ) 2 ) J 21 = (1 − m ) p ∗ 2 + β 2 (1 − m ) p ∗ 2 2 (1 + β 1 (1 − m ) n ∗ 2 + β 2 p ∗ 2 ) 2 J 22 = (1 − m ) n ∗ 2 + β 1 (1 − m ) 2 n ∗ 2 2 (1 + β 1 (1 − m ) n ∗ 2 + β 2 p ∗ 2 ) 2 − σ 1 .
Nilai eigen dari J ∗ ( ξ 6 ) mempunyai bagian real negatif jika trace ( J ∗ ( ξ 6 )) < 0 dan det ( J ∗ ( ξ 6 )) > 0 [20] . Selanjutnya, dari persamaan (15) , akan ditunjukkan bahwa:
• trace ( J ∗ ( ξ 6 )) < 0 :
trace ( J ∗ ( ξ 6 )) = J 11 + J 22 . (16)
Karena semua nilai parameter positif, maka trace ( J ∗ ( ξ 6 )) < 0 . • det ( J ∗ ( ξ 6 )) > 0 : det ( J ∗ ( ξ 6 )) = J 11 J 22 − J 21 J 12 . (17)
Karena semua nilai parameter positif, maka det ( J ∗ ( ξ 6 )) > 0 .
Dengan demikian, titik kesetimbangan ξ 6 stabil asimtotik lokal karena trace ( J ∗ ( ξ 6 )) < 0 dan det ( J ∗ ( ξ 6 )) > 0 .

## 3.3. Analisis Bifurkasi
Pada bagian ini disajikan analisis bifurkasi sistem (6) pada titik kesetimbangan ξ 4 dan ξ 6 dalam Teorema 4 beserta pembuk- tiannya. Dalam hal ini diambil parameter perlindungan mangsa m sebagai parameter bifurkasi.
Teorema 4. Sistem (6) mengalami perubahan kestabilan atau bifurkasi transcritical pada titik kesetimbangan ξ 4 ketika para- meter perlindungan mangsa m melewati titik percabangan yaitu m ∗ = 1 − σ 1 k 2 (1 − σ 1 β 1 ) .
Bukti. Berdasarkan teorema Sotomayor [21] , diambil parameter perlindungan mangsa m menjadi parameter bifurkasi dengan nilai bifurkasi m ∗ = 1 − σ 1 k 2 (1 − σ 1 β 1 ) , dengan mengevaluasi matriks Jacobian sistem (6) di ( ξ 4 , m ∗ ), yaitu:
J ( ξ 4 , m ∗ ) =   − 1 − σ 1 1 0 0 0 0 0 − σ 2   . (18)
Dari matriks Jacobian pada persamaan (18) , jelas ditunjukkan bahwa matriks J ( ξ 4 , m ∗ ) memiliki nilai eigen nol, misalkan λ 2 = 0 . Selanjutnya, perlu untuk mencari vektor eigen kanan dan kiri dari λ 2 = 0 . Misalkan u = ( u 1 , u 2 , u 3 ) T dan v = ( v 1 , v 2 , v 3 ) T masing-masing adalah vektor eigen kanan dan kiri dari λ 2 = 0 . Kedua vektor eigen tersebut dapat dicari dengan menyelesaikan J ( ξ 4 , m ∗ ) u = 0 dan J T ( ξ 4 , m ∗ ) v = 0 . Oleh karena itu, dipe- roleh u = ( σ 1 , − 1 , 0) T dan v = (0 , 1 , 0) T .
Sistem (6) dapat ditulis kembali dalam bentuk vektor pada persamaan (19) :
x ′ = f ( x ) , (19)
dimana x = ( n, p, k ) T dan
f ( x ) =     n ( 1 − n k ) − (1 − m ) np 1+ β 1 (1 − m ) n + β 2 p (1 − m ) np 1+ β 1 (1 − m ) n + β 2 p − σ 1 p σ 2 ( k − k 1 ) ( 1 − k − k 1 δ k )     . (20)
JJoM | Jambura J. Math
Untuk menggunakan teorema Sotomayor, diperlukan untuk memeriksa tiga kondisi pada persamaan (21) , yaitu:
v T f m ( ξ 4 , m ∗ ) = 0 , v T ( Df m ( ξ 4 , m ∗ ) u ) ̸ = 0 , dan
v T [ D 2 f ( ξ 4 , m ∗ ) ( u , u ) ] ̸ = 0 . (21)
• Untuk kondisi pertama, diperlukan f m sedemikian hingga diperoleh:
f m ( x ) =    (1+ β 2 p ) np (1+ β 1 (1 − m ) n + β 2 p ) 2 − (1+ β 2 p ) np (1+ β 1 (1 − m ) n + β 2 p ) 2 0
   dan
f m ( ξ 4 , m ∗ ) =   0 0 0   . (22)
oleh karena itu, dapat ditunjukkan bahwa v T f m ( ξ 4 , m ∗ ) = 0 .
• Untuk kondisi kedua, diperlukan Df m ( ξ 4 , m ∗ ) sedemikian hingga diperoleh:
Df m ( x ) =   D 11 D 12 0 D 21 D 22 0 0 0 0   (23)
dimana
D 11 = (1 + β 2 p ) p (1 + β 1 (1 − m ) n + β 2 p ) 2 − 2 β 1 (1 + β 2 p ) (1 − m ) np (1 + β 1 (1 − m ) n + β 2 p ) 3 D 12 = (1 + 2 β 2 p ) n (1 + β 1 (1 − m ) n + β 2 p ) 2 − 2 β 2 (1 + β 2 p ) np (1 + β 1 (1 − m ) n + β 2 p ) 3 D 21 = 2 β 1 (1 + β 2 p ) (1 − m ) np (1 + β 1 (1 − m ) n + β 2 p ) 3 − (1 + β 2 p ) p (1 + β 1 (1 − m ) n + β 2 p ) 2 D 22 = 2 β 2 (1 + β 2 p ) np (1 + β 1 (1 − m ) n + β 2 p ) 3 − (1 + 2 β 2 p ) n (1 + β 1 (1 − m ) n + β 2 p ) 2 .
dan
Df m ( ξ 4 , m ∗ ) =   0 k 2 (1 − σ 1 β 1 ) 2 0 0 − k 2 (1 − σ 1 β 1 ) 2 0 0 0 0   . (24)
Oleh karena itu, dapat ditunjukkan bahwa v T ( Df m ( ξ 4 , m ∗ ) u ) = k 2 (1 − σ 1 β 1 ) 2 ̸ = 0 . • Untuk kondisi ketiga, diperlukan D 2 f ( ξ 4 , m ∗ ) ( u , u ) sede- mikian hingga diperoleh:
D 2 f ( ξ 4 , m ∗ ) ( u , u ) =   − 2 σ 2 1 k 2 + 2 β 2 σ 1 (1 − σ 1 β 1 ) − 2 β 2 σ 1 (1 − σ 1 β 1 ) 0   . (25) Oleh karena itu, dapat ditunjukkan bahwa v T [ D 2 f ( ξ 4 , m ∗ ) ( u , u ) ] = − 2 β 2 σ 1 (1 − σ 1 β 1 ) ̸ = 0 .
Karena ketiga kondisi terpenuhi, maka berdasarkan teore- ma Sotomayor, titik kesetimbangan ξ 4 mengalami bifurkasi tran- scritical ketika
m ∗ = 1 − σ 1 k 2 (1 − σ 1 β 1 ) .
Ketika m > 1 −
σ 1 k 2 (1 − σ 1 β 1 ) , titik kesetimbangan ξ 4 me- rupakan titik simpul yang stabil asimtotik lokal, ketika m < 1 − σ 1 k 2 (1 − σ 1 β 1 ) , titik kesetimbangan ξ 4 menjadi titik pelana yang tidak stabil, dan pada saat bersamaan muncul ξ 6 yang stabil. Analisis bifurkasi titik kesetimbangan ξ 6 dapat dilakukan dengan menggunakan hasil analisis kestabilannya. Hal ini dapat dilakukan menggunakan matriks Jacobian pada persamaan (14) . Namun, karena komponen matriks pada persamaan (14) masih memuat n ∗ 2 dan p ∗ 2 yang ekspresinya cukup panjang dan rumit, maka analisis bifurkasi titik kesetimbangan ξ 6 terlalu rumit un- tuk dihitung secara analitik. Oleh karena itu, analisis bifurkasi titik kesetimbangan ξ 6 dilakukan secara numerik dengan meng- ontinuasi nilai parameter m .

## 3.4. Simulasi Numerik
Pada bagian ini, dilakukan simulasi numerik dari sistem per- samaan (6) untuk menggambarkan hasil analitik yang telah dila- kukan sebelumnya. Parameter yang digunakan disajikan pada Ta- bel 2 .
Tabel 2. Nilai parameter Parameter Nilai Referensi β 1 0,50 [5] β 2 0,01 [1] σ 1 1,00 [5] σ 2 0,40 [5]
Simulasi awal dilakukan menggunakan nilai parameter pa- da Tabel 2 dan nilai parameter untuk kasus daya dukung mening- kat diasumsikan yaitu m = 0 , 5 , k 1 = 8 , dan k 2 = 10 . Dite- mukan fakta bahwa sistem (6) mengalami perubahan kestabilan atau bifurkasi yang ditunjukkan pada Gambar 1 untuk kasus daya dukung meningkat. Kurva merah dan kurva biru menunjukkan keadaan tidak stabil dan stabil pada setiap titik kesetimbangan. H menunjukkan titik terjadinya bifurkasi Hopf serta BP menun- jukkan titik terjadinya bifurkasi Transcritical .
Gambar 1 menunjukkan terjadinya bifurkasi Transcritical pa- da titik m = 0 , 8 . Hal ini menunjukkan bahwa sistem mengalami perubahan kestabilan kepunahan pemangsa pada m = 0 , 8 , de- ngan kata lain, terdapat titik kesetimbangan kepunahan pemang- sa ξ 4 = ( k 2 , 0 , k 2 ) yang stabil pada interval nilai parameter m = (0 , 8; 1) dan tidak stabil pada interval m = (0; 0 , 8 ). Gambar 1 juga menunjukkan bahwa sistem mengalami bi- furkasi Hopf Supercritical pada m = 0 , 37 karena terdapat siklus
Gambar 1. Bifurkasi dengan parameter m dan kasus daya du- kung meningkat ( k 2 > k 1 )
periodik yang stabil setelah bifurkasi, dengan kata lain, terjadi ti- tik kesetimbangan hidup berdampingan ξ 6 = ( n ∗ , p ∗ , k 2 ) yang stabil pada interval nilai parameter m = (0 , 37; 0 , 8) . Adapun untuk interval nilai parameter m = (0; 0 , 37) , sistem memili- ki solusi periodik atau dengan kata lain pertumbuhan populasi mangsa dan pemangsa mengalami siklus yang berkala. Jumlah mangsa yang meningkat akan menyebabkan peningkatan jumlah pemangsa. Kemudian, dengan penurunan populasi mangsa, po- pulasi pemangsa juga dapat menurun.
Untuk memperjelas Gambar 1 , ditunjukkan simulasi da- lam potret fase mengenai bagaimana pengaruh parameter per- lindungan mangsa dengan kasus daya dukung meningkat. Si- mulasi dilakukan dengan nilai parameter pada Tabel 2 dengan 3 kasus nilai parameter perlindungan mangsa berbeda yaitu m = 0 , 1; 0 , 5; 0 , 9 serta batas dukung yang digunakan adalah k 1 = 4 , k 2 = 12 . Potret fase dengan 3 kasus nilai m beserta penjela- sannya disajikan pada Gambar 2 .
Gambar 2. Potret fase dengan m=0.1 dan kasus daya du- kung meningkat ( k 2 > k 1 )
Gambar 2 menunjukkan kasus pertama pada saat perlin- dungan mangsa kecil m = 0 , 1 dan berada pada interval m = (0; 0 , 37) , terdapat titik kesetimbangan kepunahan mangsa- pemangsa ξ 1 = (0 , 0 , k 1 ) dan ξ 2 = (0 , 0 , k 2 ) , titik kese- timbangan kepunahan pemangsa ξ 3 = ( k 1 , 0 , k 1 ) dan ξ 4 = ( k 2 , 0 , k 2 ) , serta titik kesetimbangan hidup berdampingan ξ 5 = ( n ∗ , p ∗ , k 1 ) yang tidak stabil dan merupakan titik pelana. Ke- mudian untuk titik kesetimbangan hidup berdampingan ξ 6 = ( n ∗ , p ∗ , k 2 ) tidak stabil namun terdapat limit cycle disekitar titik
ξ 6 yang artinya pertumbuhan populasi mangsa dan pemangsa te- tap ada dan mengalami siklus yang berkala, seperti ditunjukkan pada Gambar 3 .
Gambar 3. Potret fase dengan m=0.5 dan kasus daya du- kung meningkat ( k 2 > k 1 )
Gambar 3 menunjukkan kasus kedua pada saat perlindung- an mangsa m = 0 , 5 yaitu berada pada interval m = (0 , 37; 0 , 8) , terdapat titik kesetimbangan kepunahan mangsa-pemangsa ξ 1 = (0 , 0 , k 1 ) dan ξ 2 = (0 , 0 , k 2 ) , titik kesetimbangan kepunahan pemangsa ξ 3 = ( k 1 , 0 , k 1 ) dan ξ 4 = ( k 2 , 0 , k 2 ) , serta titik ke- setimbangan hidup berdampingan ξ 5 = ( n ∗ , p ∗ , k 1 ) yang tidak stabil dan merupakan titik pelana. Kemudian untuk titik kesetim- bangan hidup berdampingan ξ 6 = ( n ∗ , p ∗ , k 2 ) stabil asimtotik dan merupakan titik spiral, seperti ditunjukkan pada Gambar 4 .
Gambar 4. Potret fase dengan m=0.9 dan kasus daya du- kung meningkat ( k 2 > k 1 )
Gambar 4 menunjukkan kasus ketiga pada saat perlin- dungan mangsa besar m = 0 , 9 yaitu berada pada interval m = (0 , 8; 1) , terdapat titik kesetimbangan kepunahan mangsa- pemangsa ξ 1 = (0 , 0 , k 1 ) dan ξ 2 = (0 , 0 , k 2 ) , serta titik ke- setimbangan kepunahan pemangsa ξ 3 = ( k 1 , 0 , k 1 ) yang tidak stabil dan merupakan titik pelana. Kemudian untuk titik kesetim- bangan kepunahan pemangsa ξ 4 = ( k 2 , 0 , k 2 ) stabil asimtotik dan merupakan titik simpul. Dengan perlindungan mangsa yang besar, jelas bahwa populasi pemangsa mengalami kepunahan ka- rena populasi mangsa diberikan perlindungan yang besar sehing- ga pemangsa kesulitan untuk memangsa dan mendapat energi untuk hidup dan berkembang biak.
Simulasi selanjutnya dilakukan menggunakan nilai parame- ter pada Tabel 2 dan nilai parameter untuk kasus daya dukung menurun diasumsikan yaitu m = 0 , 2; k 1 = 8; dan k 2 = 5 . Dite- mukan fakta bahwa sistem (6) mengalami perubahan kestabilan
atau bifurkasi yang ditunjukkan pada Gambar 5 untuk kasus da- ya dukung menurun. Kurva merah dan kurva biru menunjukkan keadaan tidak stabil dan stabil pada setiap titik kesetimbangan. Dengan BP menunjukkan titik terjadinya bifurkasi Transcritical .
Gambar 5. Bifurkasi dengan parameter m dan kasus daya du- kung menurun ( k 2 < k 1 )
Gambar 5 menunjukkan terjadi bifurkasi Transcritical pada m = 0 , 6 . Hal ini menunjukkan bahwa sistem mengalami peru- bahan kestabilan kepunahan pemangsa pada m = 0 , 6 , atau de- ngan kata lain, terdapat titik kesetimbangan kepunahan pemang- sa ξ 4 = ( k 2 , 0 , k 2 ) yang stabil pada interval nilai parameter m = (0 , 6; 1) dan tidak stabil ketika m = (0; 0 , 6) . Kemudian terda- pat titik kesetimbangan hidup berdampingan ξ 6 = ( n ∗ , p ∗ , k 2 ) yang stabil pada interval nilai parameter m = (0; 0 , 6) .
Simulasi mengenai pengaruh parameter perlindungan mangsa dengan kasus daya dukung menurun dapat dilakukan de- ngan membagi menjadi 2 kasus nilai perlindungan mangsa m , ya- itu pada interval m = (0; 0 , 6) dan m = (0 , 6; 1) . Secara ku- alitatif, penjelasan mengenai pengaruh parameter perlindungan mangsa dengan kasus daya dukung menurun sama dengan kasus daya dukung meningkat, yang ditampilkan pada Gambar 3 dan Gambar 4 .

## 4. Kesimpulan
Pada penelitian ini telah dilakukan pencarian titik kesetim- bangan, analisis kestabilan, dan simulasi numerik sehingga di- peroleh bahwa parameter perlindungan mangsa dan daya du- kung variabel berpengaruh terhadap kestabilan sistem mangsa- pemangsa dengan fungsi respon Beddington-DeAngelis serta dite- mukan adanya bifurkasi atau perubahan kestabilan pada nilai pa- rameter perlindungan mangsa tertentu pada saat dilakukan si- mulasi numerik. Untuk kasus daya dukung meningkat ( k 2 > k 1 ) , terdapat titik kesetimbangan kepunahan pemangsa ξ 4 = ( k 2 , 0 , k 2 ) yang stabil untuk nilai parameter perlindungan mang- sa yang besar yaitu pada rentang (0,8; 1). Terdapat titik kesetim- bangan hidup berdampingan ξ 6 = ( n ∗ , p ∗ , k 2 ) yang stabil untuk nilai parameter perlindungan mangsa pada rentang (0,37; 0,8). Untuk nilai parameter perlindungan mangsa yang kecil yaitu pada rentang (0; 0,37), populasi mangsa dan pemangsa bersifat perio- dik atau terjadi fluktuasi dengan kata lain pertumbuhan popula- si mangsa dan pemangsa mengalami siklus yang berkala. Untuk
kasus daya dukung menurun atau ( k 2 < k 1 ) , sistem mencapai kesetimbangan kepunahan pemangsa yang stabil untuk nilai pa- rameter perlindungan mangsa yang besar. Sistem juga mencapai kesetimbangan hidup berdampingan yang stabil untuk nilai pa- rameter perlindungan mangsa pada rentang (0; 0,6). Perlindung- an mangsa memberikan dampak pada upaya konservasi populasi mangsa dan pemangsa. Ketika perlindungan mangsa kecil, akan terjadi fluktuasi atau naik turunnya pertumbuhan populasi mang- sa dan pemangsa, keadaan tersebut lebih rentan oleh pengaruh eksternal seperti cuaca, jika cuaca tidak mendukung dapat meng- akibatkan terjadinya kepunahan. Dengan adanya perlindungan mangsa yang lebih besar sistem dapat lebih stabil dan tidak mu- dah terganggu oleh faktor-faktor eksternal yang dapat menye- babkan kepunahan.
Kontribusi Penulis. Aisyiah Kholifatul Hanifah: Konseptualisasi, pe- ngumpulan data, metodologi, analisis formal, investigasi, penulisan draft, dan visualisasi. Abadi: Konseptualisasi, validasi, penulisan pem- bahasan dan pengeditan, supervisi. Semua penulis telah membaca dan menyetujui versi manuskrip yang diterbitkan.
Ucapan Terima Kasih. Para penulis menyampaikan terima kasih kepada editor dan reviewer atas pembacaan yang cermat, kritik yang mendalam, dan rekomendasi yang praktis untuk meningkatkan kualitas tulisan ini.
Pembiayaan. Penelitian ini tidak menerima pembiayaan eksternal.
Konflik Kepentingan. Para penulis menyatakan tidak ada konflik kepen- tingan yang terkait dengan artikel ini.

## Referensi
[1] J. P. Tripathi, S. Abbas, and M. Thakur, “Dynamical analysis of a prey– predator model with Beddington–DeAngelis type function response incor- porating a prey refuge,” Nonlinear Dyn. , vol. 80, no. 1–2, pp. 177–196, 2014, doi: 10.1007/s11071-014-1859-2 .
[2] S. Chakraborty, P. K. Tiwari, S. K. Sasmal, S. Biswas, S. Bhattacharya, and J. Chattopadhyay, “Interactive effects of prey refuge and additional food for predator in a diffusive predator-prey system,” Appl. Math. Model. , vol. 47, pp. 128–140, 2017, doi: 10.1016/j.apm.2017.03.028 .
[3] S. Samanta, R. Dhar, I. M. Elmojtaba, and J. Chattopadhyay, “The Role of Additional Food in A Predator-prey Model with A Pey Refuge,” J. Biol. Syst. , vol. 24, no. 2–3, pp. 345–365, 2016, doi: 10.1142/S0218339016500182 .
[4] J. H. Connell, S. E. Monographs, N. Winter, and J. H. Connell, “A Predator- Prey System in the Marine Intertidal Region . I . Balanus glandula and Several Predatory Species of Thais,” Ecol. Soc. Am. , vol. 40, no. 1, pp. 49–78, 1970.
[5] N. Al-Salti, F. Al-Musalhi, V. Gandhi, M. Al-Moqbali, and I. Elmojtaba, “Dyna- mical analysis of a prey-predator model incorporating a prey refuge with variable carrying capacity,” Ecol. Complex. , vol. 45, no. July 2020, p. 100888, 2021, doi: 10.1016/j.ecocom.2020.100888 .
[6] E. González-Olivares and R. Ramos-Jiliberto, “Dynamic consequences of prey refuges in a simple model system: More prey, fewer predators and enhan- ced stability,” Ecol. Modell. , vol. 166, no. 1–2, pp. 135–146, 2003, doi: 10.1016/S0304-3800(03)00131-5 .
[7] Y. Huang, F. Chen, and L. Zhong, “Stability analysis of a prey-predator model with holling type III response function incorporating a prey re- fuge,” Appl. Math. Comput. , vol. 182, no. 1, pp. 672–683, 2006, doi: 10.1016/j.amc.2006.04.030 .
[8] Z. Ma, S. Wang, W. Li, and Z. Li, “The effect of prey refuge in a patchy predator-prey system,” Math. Biosci. , vol. 243, no. 1, pp. 126–130, 2013, doi: 10.1016/j.mbs.2013.02.011 .
[9] D. Mukherjee, “The effect of refuge and immigration in a predator-prey sys- tem in the presence of a competitor for the prey,” Nonlinear Anal. Real World Appl. , vol. 31, pp. 277–287, 2016, doi: 10.1016/j.nonrwa.2016.02.004 .
[10] A. Sih, “Prey refuges and predator-prey stability,” Theor. Popul. Biol. , vol. 31, no. 1, pp. 1–12, 1987, doi: 10.1016/0040-5809(87)90019-0 .
[11] N. M. S. Sundari and M. Valliathal, “Analysing a prey predator model for
stability with prey refuge-stage structure on predator,” J. Phys. Conf. Ser. , vol. 1139, no. 1, 2018, doi: 10.1088/1742-6596/1139/1/012016 .
[12] M. Moustafa, M. H. Mohd, A. I. Ismail, and F. A. Abdullah, “Dynamical analysis of a fractional-order Rosenzweig–MacArthur model incorporating a prey refuge,” Chaos, Solitons and Fractals , vol. 109, pp. 1–13, 2018, doi:
10.1016/j.chaos.2018.02.008 .
[13] M. K. A. Al-Moqbali, N. S. Al-Salti, and I. M. Elmojtaba, “Prey-predator models with variable carrying capacity,” Mathematics , vol. 6, no. 6, Jun. 2018, doi: 10.3390/math6060102 .
[14] C. Ganguli, T. K. Kar, and P. K. Mondal, “Optimal harvesting of a prey-predator model with variable carrying capacity,” Int. J. Biomath. , vol. 10, no. 5, pp. 1– 24, 2017, doi: 10.1142/S1793524517500693 .
[15] P. S. Meyer and J. H. Ausubel, “Carrying capacity: A model with logistically varying limits,” Technol. Forecast. Soc. Change , vol. 61, no. 3, pp. 209–214, 1999, doi: 10.1016/S0040-1625(99)00022-0 .
[16] C. S. Holling, “The Components of Predation as Revealed by a Study of Small- Mammal Predation of the European Pine Sawfly,” Can. Entomol. , vol. 91, no. 547, pp. 293–320, 1959, doi: 10.2351/1.5058368 .
[17] M. P. Hassell, “Mutual Interference between Searching Insect Parasites,” J. Anim. Ecol. , vol. 40, no. 2, p. 473, 1971, doi: 10.2307/3256 .
[18] J. R. Beddington, “Mutual Interference Between Parasites or Predators and its Effect on Searching Efficiency,” J. Anim. Ecol. , vol. 44, no. 1, p. 331, 1975, doi: 10.2307/3866 .
[19] D. L. DeAngelis, R. A. Goldstein, and R. V. O’Neill, “A Model for Tropic Inte- raction,” Ecology , vol. 56, no. 4, pp. 881–892, 1975, doi: 10.2307/1936298 .
[20] W. E. Boyce, R. C. DiPrima, and D. B. Meade, Elementary differential equations and boundary value problems , 12th Edition. John Wiley & Sons, 2021.
[21] L. Perko, Equations and Dynamical Systems . New York: Springer, 2001, doi: 10.1007/978-1-4613-0003-8 .