Limits: Journal of Mathematics and Its Applications E-ISSN: 2579-8936 P-ISSN: 1829-605X Vol. 17, No. 1, Juli 2020, 83-95 DOI: http://dx.doi.org/10.12962/limits.v17i1.6753

## Model Matematika Pemanenan Ikan dengan Kebijakan Panen Selektif
Lisa Risfana Sari 1 *
1 Universitas Internasional Semen Indonesia, Gresik Indonesia e-mail: lisa.sari@uisi.ac.id
Diajukan:9 April 2020,  Diperbaiki: 17 Juni 2020, Diterima:21 Juli 2020

## Abstrak
Pengelolaan  sumber  daya  perikanan  termasuk  sebagai  aspek  penting  yang  dipertimbangkan  oleh negara maritim. Kelimpahan sumber daya ikan dapat dipertahankan dengan strategi penangkapan ikan yang tepat, salah satunya adalah kebijakan panen selektif. Dalam studi ini, dinamika kepadatan populasi ikan dipelajari  menggunakan  model predator-prey yang  dimodifikasi.  Proses  panen  selektif  yang memperhitungkan  usia atau ukuran ikan siap panen dinyatakan  sebagai  penundaan  waktu  dalam  proses panen.  Analisis  model  dilakukan  dengan  menentukan  titik  kesetimbangan  model  dan  stabilitas  titik keseimbangan model. Ada empat titik keseimbangan model, yang mewakili kondisi trivial, kepunahan prey, kepunahan predator , dan eksistensi predator-prey . Perilaku dinamis model diilustrasikan melalui simulasi numerik  dengan  beberapa  skenario.  Hasil  simulasi  numerik  menunjukkan  bahwa  waktu  tunda mempengaruhi stabilitas beberapa titik setimbang, sehingga menghasilkan dinamika populasi yang lebih beragam.
Kata Kunci : pemanenan ikan, p redator-prey , simulasi numerik, waktu tunda.

## Abstract
The abundance of fish resources can be maintained with the right fishing strategies, one of which is the selective harvest policy. In this study, the dynamics of fish population density were studied using a modified predator-prey model. The selective harvest process that takes into account the age or size at which fish are ready to harvest is expressed as a time delay in the harvesting process. The model analysis involves determining the model's equilibrium point and the model's equilibrium point stability. There are four equilibrium points of the model, which represent trivial conditions, prey extinction, predator extinction, and predator-prey existence. The dynamic behavior of the model is illustrated through numerical simulations with several scenarios. Numerical simulation results show that time delay affects the stability of several equilibrium points, resulting in a more diverse population dynamics.
Keywords: fish harvesting, predator-prey, numerical simulation, time delay.

## 1 Pendahuluan
Ekologi  akuatik  adalah  sub-disiplin  ekologi  yang  mempelajari  interaksi  atau  keterkaitan organisme  yang  hidup  di  perairan  dengan  lingkungannya.  Termasuk  didalamnya  ialah  lautan, muara,  danau,  kolam,  dan  sungai.  Lingkungan  alami  mempengaruhi  pertumbuhan  dan
perkembangan  organisme  hidup.  Tentu  saja,  hal  ini  berkontribusi  pada  kelangsungan  hidup organismenya [1].
Ikan adalah salah satu organisme akuatik. Ekologi akuatik memiliki posisi yang berpengaruh di  bidang  perikanan.  Manfaat  mempelajari  ekologi  perairan  adalah  mengetahui  berbagai  jenis ekosistem perairan, juga parameter atau faktor pembatas lingkungan bagi organisme yang hidup di dalamnya [1]. Hal ini dapat digunakan untuk mendukung manajemen sumber daya perikanan.
Sumber daya perikanan berpotensi untuk meningkatkan pertumbuhan dan kekayaan negara. Penerapan  manajemen  perikanan  yang  efektif  dapat  meningkatkan  produktivitas  perikanan. Kelimpahan  sumber  daya  perikanan  yang  dikelola  dengan  baik  dapat  menjamin  pemanfaatan berkelanjutan untuk kesejahteraan masyarakat. Oleh karena itu, diperlukan kebijakan manajemen sumber daya untuk menghindari eksploitasi ikan dan menjaga cadangan ikan.
Pemanenan  ikan  mempengaruhi  kelimpahan  sumber  daya  ikan  dan  hasil  perikanan maksimum. Kebijakan pemanenan selektif adalah upaya untuk melestarikan sumber daya ikan. Dalam pemanenan selektif, ikan yang siap panen adalah ikan yang telah mencapai ukuran atau usia standar tertentu [2]–[5]. Hal ini dilakukan untuk memastikan keberadaan, ketersediaan, dan keberlanjutan sumber daya ikan sekaligus mempertahankan dan meningkatkan kualitas nilai dan keanekaragaman.
Pemodelan matematika menggunakan model predator-prey telah banyak digunakan dalam studi ekologi akuatik. Model predator-prey digunakan untuk mempelajari interaksi dua spesies ikan dalam suatu ekosistem yang dibagi menjadi spesies mangsa ( prey ) dan spesies predator . Hasil beberapa penelitian menunjukkan bahwa model dapat menggambarkan dinamika populasi ikan, menjelaskan pengaruh berbagai faktor pada populasi, serta memprediksi perilaku dinamis populasi [5]–[9].
Berkaitan dengan manajemen perikanan, penelitian ini dilakukan untuk menguji bagaimana kebijakan  panen  selektif  mempengaruhi  kelimpahan  sumber  daya  ikan  dari  perspektif  model matematika.  Dalam  studi  ini,  model predator-prey yang  dimodifikasi  diformulasikan  dengan fungsi  respon  Beddington-DeAngelis  dalam  pemangsaan  dan  penundaan  panen  pada  mangsa ( prey ), yang memperhitungkan pemanenan ikan secara selektif di atas usia atau ukuran tertentu.
Analisis dinamis model dilakukan untuk melacak kondisi keseimbangan dalam jangka panjang. Pada  bagian  terakhir,  pengaruh  seleksi  pada  panen  terhadap  hasil  perikanan  dibahas  dalam serangkaian simulasi.

## 2 Model Matematika
Pemodelan  matematika  dikembangkan  berdasarkan  model predator-prey yang  diusulkan oleh [5], [6], [10]–[14]. Model dinyatakan dalam sistem diferensial tundaan nonlinier dengan dua variabel  yang  merupakan  fungsi  terhadap  waktu.  Fungsi  tersebut  mewakili  populasi  kepadatan mangsa  ( prey )  dan  populasi  kepadatan predator .  Model predator-prey yang  dimodifikasi dijelaskan dalam sistem pada Persamaan (1).
      2 1 1 2 2 ( ) ( ) ( ) ( ), ( ) ( ) ,
( ) ( ) ( ) ( ) . ( ) dx t r x t px t f x t y t H x t dt dy t y t y t r H y t dt x t k 
               
(1)
Dalam sistem (1), parameter 1 r dan 2 r masing-masing menunjukkan tingkat pertumbuhan intrinsik mangsa dan predator . Parameter p menunjukkan kekuatan persaingan di antara individu mangsa,  menunjukkan  nilai  maksimum  tingkat  pengurangan  per  kapita predator ,  dan k mengukur tingkat perlindungan lingkungan terhadap predator .
Pemangsaan prey oleh predator digambarkan sebagai fungsi respon Beddington-DeAngelis pada  Persamaan  (2).  Pemilihan  fungsi  respon  ini  dilatarbelakangi  asumsi  bahwa  populasi prey tersaturasi, dimana laju pemangsaan akan menurun seiring dengan bertambahnya populasi baik prey maupun predator . Hal ini sesuai dengan model yang diusulkan oleh [11], [12], [14].
  ( ) ( ) ( ), ( ) ( ) ( ) x t y t f x t y t a bx t cy t     . (2)
Parameter  adalah  nilai  maksimum  yang  tingkat  pengurangan prey per  kapita  dan a merupakan  tingkat  perlindungan  lingkungan  untuk  mempertahankan  kelangsungan  hidup prey . Sementara  itu,  parameter b dan c koefisien  terkait  dengan  kontak  mutual  diantara prey dan predator .
Populasi  ikan  termasuk prey dan predator dipanen  menggunakan  fungsi  panen  linier. Diasumsikan predator adalah ikan yang berukuran lebih besar atau berusia lebih dewasa daripada spesies prey .  Pemanenan  selektif  memastikan  pembatasan  untuk  menangkap  ikan  sebelum mencapai ukuran atau usia tertentu [3]. Ini diungkapkan dengan memasukkan penundaan waktu dalam  masa  panen  [2],  [4]–[6].  Didasari  oleh  hal  tersebut,  kebijakan  pemanenan  selektif  pada model hanya diterapkan untuk spesies prey . Fungsi panen untuk kedua spesies dijelaskan sebagai berikut.
1 1 1 2 2 2 ( ( )) ( ) ( ( )) ( )
H x t q E x t
H y t q E y t       (3)
Parameter 1 q and 2 q masing-masing  menunjukkan  koefisien  penangkapan  dari prey dan predator . Sementara itu, parameter 1 E dan 2 E menunjukkan upaya panen untuk populasi prey dan predator .

## 3 Hasil dan Pembahasan
Pada bagian ini akan dibahas hasil analisis numerik beserta simulasinya. Kedua hasil tersebut digunakan untuk menggambarkan perilaku dinamis populasi.
3.1 Titik Kesetimbangan
Titik kesetimbangan model diperoleh dari solusi setimbang saat ( ) 0 dx t dt  dan ( ) 0 dy t dt  . Diperoleh empat titik kesetimbangan model yaitu;
1. titik trivial 0 (0,0) I  ,
2. titik kepunahan prey
  2 2 2 1 0, k r q E I      
   
 ,
3. titik kepunahan predator 1 1 1 2 , 0 r q E I p         ,
4. titik eksistensi prey dan predator    2 2 2 3 , r q E k I               ,
dimana  adalah akar positif riil dari persamaan berikut
     
2 2 2 1 2 1 2 1 2 1 2 P p c b ap b ckp c ck a k                             
dengan
1 1 1 1 2 2 2 2 r q E r q E      
Titik kesetimbangan mewakili kepadatan populasi prey dan predator pada masa setimbang. Titik
0 I menggambarkan kepunahan bagi prey dan predator . Titik ini selalu eksis. Titik 1 I dan 2 I menggambarkan  kepunahan  salah  satu  spesies  sementara  yang  lain  bertahan.  Titik 1 I dan 2 I masing-masing  eksis  jika  hanya  jika 2 2 2 0 r q E   dan 1 1 1 0 r q E   .  Sementara  itu,  titik 3 I menjelaskan  bahwa  kedua  spesies  ini  bertahan  hidup.  Titik  ini  eksis  jika  dan  hanya  jika
2 2 2 0 r q E   dan 1 1 1 0 r q E   .

## 3.2 Kestabilan Lokal Titik Kesetimbangan
Linearisasi sistem persamaan diferensial tundaan diperlukan dalam analisis kestabilan titik kesetimbangan.  Persamaan  karakteristik  dari  sistem  terlinearisasi  model  (1)  dinyatakan  dalam Persamaan (4).
11 12
21 22 0 j j j j      (4)
dimana
        * * * * 11 1 1 1 * * 2 * * * * * 12 * * 2 * * 2 * 21 2 * * 22 2 2 2 * 2 2 y bx y j r px e q E a bx cy a bx cy x cx y j a bx cy a bx cy y j
x k y j r q E x k     
 
                
 
   
Penentuan kestabilan lokal pada model sistem persamaan diferensial tundaan tidak mudah. Oleh  karenanya,  digunakan  analisis  numerik  sebagai  acuan.  Penentuan  nilai  parameter  pada analisis tersebut merujuk pada [5]–[9], [15].
Berkaitan  dengan  syarat  eksistensi  titik  kesetimbangan,  diperoleh  empat  macam  kondisi yang  digunakan  untuk  menentukan  dan  mengamati  kestabilan  titik  kesetimbangan.  Keempat kondisi tersebut dijelaskan berikut.
1. Kondisi 1 1 1 r q E  dan 2 2 2 r q E 
Kondisi ini  mewakili keadaan over  eksploitasi  kedua  spesies,  dengan laju pertumbuhan lebih  kecil  dibandingkan  upaya  penangkapan  ikan.  Pada  kondisi  ini  digunakan  nilai parameter 3 1 5 10 , r    3 1 0.1 10 , q    1 72, E  3 2 8 10 , r    3 2 0.5 10 , q    dan 2 22.
E  Terdapat satu titik kesetimbangan eksis yaitu 0 I . Persamaan karakteristik pada titik 0 (0,0)
I  adalah

   2 3 6 3 6 ( ) 2 10 15 10 7.2 10 21.6 10 0 p e                     . (5) Berdasarkan Persamaan (5) untuk 0   , diperoleh bagian riil akar persamaan karakteristik bernilai  negatif,  sehingga  titik 0 I bersifat  stabil  asimtotik  lokal.  Sementara  itu,  untuk
0   , diperoleh 0 155.0260787   . Nilai kritis ini merupakan batas perubahan kestabilan. Kestabilan  titik  kesetimbangan  sistem  dijamin  jika ( ) 0 Re   untuk  semua  akar karakteristik.  Titik  kesetimbangan  kehilangan  stabilitasnya  hanya  jika  nilai  melewati nilai kritis tundaan yang diperoleh dari  ( ) 0 p i   . Perubahan kestabilan bisa terjadi antara kondisi stabil dan tidak stabil atau sebaliknya [16]. Perubahan kestabilan ini lebih lanjut akan digambarkan melalui simulasi numerik.
2. Kondisi 1 1 1 r q E  dan 2 2 2 r q E 
Situasi yang terjadi pada kondisi berikut adalah over eksploitasi hanya terjadi pada spesies prey . Terdapat dua titik kesetimbangan yang eksis, yaitu 0 I dan 1 I . Disini digunakan nilai
parameter  yaitu 3 1 5 10 , r    3 1 0.3 10 , q    1 22, E  3 2 8 10 , r    3 2 0.5 10 , q    dan
2 12. E  Substitusi masing-masing titik kesetimbangan pada persamaan (4) menghasilkan persamaan karakteristik berikut
    2 3 5 3 6 7 10 10 6.6 10 13.2 10 0 e                  untuk 0 I (6)  
2 6 6 0.0029 9.9 10 (0.0066 13.2 10 ) 0 e               untuk 1 I . (7)
Solusi Persamaan (6) dan (7) untuk 0   menunjukkan bahwa titik 0 I tidak stabil, dimana terdapat akar dengan bagian riil positif, sedangkan 1 I stabil asimtotik lokal. Lebih lanjut untuk 0   ,  hasil  menunjukkan  tidak  terdapat  nilai  kritis  tundaan.  Akibatnya,  tidak ditemukan terjadinya perubahan kestabilan pada kedua titik kesetimbangan.
3. Kondisi 1 1 1 r q E  dan 2 2 2 r q E 
Berbeda  dengan  kondisi  kedua,  pada  kondisi  ketiga  terjadi  over  eksploitasi  spesies predator .  Titik  kesetimbangan 0 I dan 2 I eksis  pada  kondisi  tersebut.  Pengamatan
dilakukan  dengan  nilai  parameter 3 1 6 10 , r    3 1 0.3 10 , q    1 12, E  3 2 3 10 , r    3 2 0.2 10 , q    dan 2 22.
E  Persamaan  karakteristik  untuk  masing-masing  titik kesetimbangan adalah     2 3 6 3 6 4.6 10 8.4 10 5.04 10 3.6 10 0 e                   untuk 0 I (8)   2 6 6 0.002 1.7 10 (0.0036 5 10 ) 0 e               untuk 2 I . (9)
Perhitungan dengan 0   untuk Persamaan (8) dan (9) mengungkapkan bahwa 0 I tidak stabil dan 2 I stabil asimtotik lokal. Kemudian untuk 0   , solusi Persamaan (8) dan (9)
mengkonfirmasi bahwa tidak terdapat nilai kritis tundaan. Dengan kata lain, tidak terjadi perubahan kestabilan pada kedua titik kesetimbangan.
4. Kondisi 1 1 1 r q E  dan 2 2 2 r q E 
Kondisi  terakhir  ini  menggambarkan  laju  pertumbuhan  kedua  spesies  ikan  masih  lebih tinggi  dibandingkan  upaya  penangkapannya.  Hal  ini  juga  mewakili  kondisi  ketika kelimpahan ikan untuk kedua spesies dapat dipertahankan. Keempat titik kesetimbangan eksis  pada  kondisi  tersebut.  Nilai  parameter  yang  digunakan  untuk  analisis  adalah
3
1 6 10 , r    3 1 0.3 10 , q    1 12, E  3 2 7 10 , r    3 2 0.2 10 , q    dan 2 22. E  Untuk masing-masing titik, persamaan karakteristiknya adalah
    2 3 6 3 6 8.6 10 15.6 10 9.36 10 3.6 10 0 e                   untuk 0 I (10)     2 3 6 3 6 3.6 10 3.2 10 3.6 10 9.4 10 0 e                   untuk 1 I (11)     2 3 6 3 6 3.4 10 15.6 10 3.4 10 9.4 10 0 e                   untuk 2 I (12)     2 3 6 3 6 3.6 10 3.2 10 3.6 10 9.4 10 0 e                   untuk 3 I (13)
Perolehan  solusi  untuk 0   pada  Persamaan  (10)-(13)  menunjukkan  bahwa  titik kesetimbangan 0 I , 1 I ,  dan 2 I tidak  stabil.  Sementara  itu,  titik 3 I stabil  asimtotik  lokal.
Analisis  dengan 0   mengungkapkan  tidak  terdapat  nilai  kritis  tundaan  sehingga perubahan kestabilan tidak terjadi.
3.3 Simulasi Numerik
Simulasi dilakukan pada keempat kondisi, dengan bantuan software Matlab. Studi numerik menunjukkan  bahwa  eksploitasi  berlebihan  dapat  menyebabkan  kepunahan  populasi  ikan. Diilustrasikan  pada  Gambar  1,  jika 1 1 1 r q E  dan 2 2 2 r q E  maka  semua  solusi  Sistem  (1) konvergen ke 0 I untuk 0   .
Gambar 1. Potret fase untuk 1 1 1 r q E  dan 2 2 2 r q E  saat 0  
Untuk  mencegah  situasi  tersebut,  dibutuhkan  penundaan  waktu  panen  pada prey yang cukup  besar.  Tampak  bahwa  terdapat  perubahan  kestabilan 0 I dengan  peningkatan  (Lihat Gambar 2). Dalam hal ini, untuk nilai tunda  yang lebih besar dari nilai kritis 0  , populasi prey dapat dilindungi dari kepunahan.
(i) 144   (ii) 154   (iii) 155   (iv) 156   Gambar 2. Solusi sistem untuk 1 1 1 r q E  dan 2 2 2 r q E 
Gambar 3 menunjukkan perilaku osilasi pada kepadatan populasi  mangsa (prey) . Osilasi semakin  besar  ketika  nilai  lebih  besar  dari  155.  Hal  ini  mengindikasikan  bahwa 155   merupakan  suatu  nilai  kritis  tundaan  yang  mewakili  waktu  tunda  panen  minimum  untuk menghambat laju penurunan populasi ikan. Waktu tunda panen dapat dikaitkan dengan kebijakan panen  selektif  berdasarkan  ukuran  atau  usia  siap  panen  ikan.  Namun  demikian,  osilasi  terjadi diantara nilai nol sehingga hasil ini kurang dapat mewakili kondisi nyata kepadatan populasi ikan.
(i) 144   (ii) 154  
(i) 155   (ii) 156   Gambar 3. Kepadatan populasi ikan untuk 1 1 1 r q E  dan 2 2 2 r q E  Gambar 4 mengungkapkan jika 1 1 1 r q E  dan 2 2 2 r q E  maka solusi Sistem (1) konvergen ke 1 I untuk 0   . Hal ini didukung Gambar 5 yang menunjukkan bahwa walaupun terjadi osilasi pada kepadatan populasi mangsa saat 0   , besar osilasi semakin kecil dan menjadi setimbang pada nilai nol. Sementara itu, jika 2 2 2 r q E  dan 1 1 1 r q E  maka solusi Sistem (1) konvergen ke
2 I untuk 0   (lihat Gambar 6). Ketiganya mewakili kondisi kepunahan salah satu spesies ikan.
Gambar 4. Solusi sistem untuk 1 1 1 r q E  dan 2 2 2 r q E  untuk  
0,60,100,144,160  
(i) 60   (ii) 100  
(i) 144   (ii) 160   Gambar 5. Kepadatan populasi ikan untuk 1 1 1 r q E  dan 2 2 2 r q E 
Gambar 6. Solusi sistem untuk 1 1 1 r q E  dan 2 2 2 r q E  untuk  

## 0,60,100,144,160
 
Eksploitasi berlebihan terhadap suatu spesies memicu penurunan kepadatan populasi spesies tersebut  dan  bahkan  dapat  menyebabkan  kepunahan.  Tampak  pada  Gambar  4  dan  6  bahwa penundaan  panen  tidak  mempengaruhi  perubahan  stabilitas  sehingga  adanya  kebijakan  panen selektif pun tidak mampu menghambat kepunahan spesies tersebut.
Kondisi dengan populasi ikan tidak tereksploitasi yaitu 1 1 1 r q E  dan 2 2 2 r q E  ditunjukkan
pada  Gambar  7.  Diilustrasikan  bahwa  solusi  Sistem  (1)  konvergen  ke 3 I untuk 0   .  Dalam
situasi ini, kedua spesies dapat dilindungi dari kepunahan. Hal ini dikarenakan populasi keduanya tetap ada atau mampu bertahan hingga pada kondisi setimbang.
Gambar 7. Solusi sistem untuk 1 1 1 r q E  dan 2 2 2 r q E  untuk   0,60,100,144,160  
Eksploitasi  berlebihan  dalam  sumber  daya  perikanan  terjadi  ketika  sumber  daya  dipanen pada  tingkat  yang  lebih  cepat  daripada  yang  laju  regenerasi.  Eksploitasi  berlebihan  merusak ekosistem dan dapat menyebabkan degradasi. Degradasi terjadi ketika net supply ekosistem rusak sehingga tidak dapat pulih sendiri dalam periode yang wajar setelah tindakan perusakan dihentikan [17].  Untuk  menjaga  keanekaragaman  hayati,  eksploitasi  berlebihan  harus  dikontrol  oleh kebijakan pemanenan yang tepat.

## 4 Simpulan
Dalam penelitian ini, dampak kebijakan panen selektif terhadap kelimpahan sumber daya ikan  diselidiki.  Model  dirumuskan  untuk  memahami  mekanisme  panen  ikan  yang  tepat  untuk menjaga stok ikan. Mekanisme tersebut terkait dengan tindakan pencegahan eksploitasi sumber daya ikan yang berlebihan.
Eksploitasi  berlebihan  untuk  kedua  spesies  membahayakan  ekosistem  dan  mengancam kelimpahan  sumber  daya  perikanan.  Membatasi  jumlah  penangkapan  ikan  atau  panen  dengan upaya pengelolaan tangkapan adalah salah satu cara untuk melindungi stok ikan dari eksploitasi berlebihan. Untuk menjaga stok ikan, waktu tunda panen dengan durasi waktu tertentu diperlukan. Hal ini diwakili oleh nilai tunda kritis. Namun, ketika eksploitasi berlebihan hanya terjadi pada salah satu spesies, penundaan panen tidak dapat mencegah penurunan kepadatan populasi sehingga stok  ikan  masih  berisiko  habis.  Oleh  karena  itu,  untuk  mengatasi  eksploitasi  sumber  daya perikanan  yang  berlebihan,  kebijakan  pemanenan  selektif  juga  perlu  didukung  oleh  kebijakan upaya penangkapan ikan terbatas.
Meskipun  model  tersebut  telah  memberikan  interpretasi  yang  jelas  tentang  berbagai fenomena  yang  terkait  dengan  pengelolaan  sumber  daya  perikanan,  model  ini  mungkin  belum memenuhi  batasan  positif  untuk  solusi.  Analisis  stabilitas  global  dari  keseimbangan  belum dipelajari  dalam  analisis  model.  Selain  itu,  waktu  tunda  panen  untuk predator perlu dipertimbangkan dalam model.

## 5 Ucapan Terima Kasih
Ucapan terima kasih diberikan kepada Universitas Internasional Semen Indonesia yang telah mendukung penelitian ini melalui Hibah Riset Bersaing UISI 2019.

## 6 Daftar Pustaka
[1] H. Latuconsina and U. G. M. Press, Ekologi Perairan Tropis: Prinsip Dasar Pengelolaan Sumber Daya Hayati Perairan . UGM PRESS, 2019.
[2] A. Mukhopadhyay, J. Chattopadhyay, and P. K. Tapaswi, “Selective harvesting in a two species fishery model,” Ecol. Modell. , 1997, doi: 10.1016/S0304-3800(96)00019-1.
[3] M. Breen, N. Graham, M. Pol, P. He, D. Reid, and P. Suuronen, “Selective fishing and balanced harvesting,” Fish. Res. , 2016, doi: 10.1016/j.fishres.2016.03.014.
[4] D. Pal, G. P. Samanta, and G. S. Mahapatra, “Selective harvesting of two competing fish species in the presence of toxicity with time delay,” Appl. Math. Comput. , 2017, doi: 10.1016/j.amc.2017.05.069.
[5] Y. Pei, M. Chen, X. Liang, and C. Li, “Model-based on fishery management systems with selective harvest policies,” Math. Comput. Simul. , 2019, doi:
10.1016/j.matcom.2018.08.009.
[6] T. K. Kar, “Selective harvesting in a prey-predator fishery with time delay,” Math. Comput. Model. , 2003, doi: 10.1016/S0895-7177(03)90099-9.
[7] Y. Lv, R. Yuan, and Y. Pei, “A prey-predator model with harvesting for fishery resource with reserve area,” Appl. Math. Model. , 2013, doi: 10.1016/j.apm.2012.07.030.
[8] M. Li, B. Chen, and H. Ye, “A bioeconomic differential algebraic predator–prey model with nonlinear prey harvesting,” Appl. Math. Model. , 2017, doi:
10.1016/j.apm.2016.09.029.
[9] S. G. Mortoja and P. Panja, “Dynamics of a predator-prey model with stage-structure on both species and anti-predator behavior,” Informatics Med. Unlocked , vol. 10, pp. 50–57, 2018, doi: 10.1016/j.imu.2017.12.004.
[10]  P. H. LESLIE, “SOME FURTHER NOTES ON THE USE OF MATRICES IN POPULATION MATHEMATICS,” Biometrika , 1948, doi: 10.1093/biomet/35.3-4.213.
[11]  J. R. Beddington, “Mutual Interference Between Parasites or Predators and its Effect on Searching Efficiency,” J. Anim. Ecol. , 1975, doi: 10.2307/3866.
[12]  D. L. DeAngelis, R. A. Goldstein, and R. V. O’Neill, “A Model for Tropic Interaction,” Ecology , 1975, doi: 10.2307/1936298.
[13]  M. A. Aziz-Alaoui and M. Daher Okiye, “Boundedness and Global Stability for a Predator-Prey Model with Modified Leslie-Gower and Holling-Type II Schemes,” Appl. Math. Lett. , 2003, doi: 10.1016/S0893-9659(03)90096-6.
[14]  S. Yu, “Global stability of a modified Leslie-Gower model with Beddington-DeAngelis functional response,” Adv. Differ. Equations , vol. 2014, no. 1, pp. 1–14, 2014, doi: 10.1186/1687-1847-2014-84.
[15]  F. Fitriah, A. Suryanto, and N. Hidayat, “Numerical Study of Predator-Prey Model with Beddington-DeAngelis Functional Response and Prey Harvesting,” J. Trop. Life Sci. , vol.
5, no. 2, pp. 105–109, May 2015, doi: 10.11594/jtls.05.02.09.
[16]  Q. Gao and H. R. Karimi, Stability, control and application of time-delay systems . Elsevier, 2019.
[17]  K. Steier, Gabriela, Patel, Ed., International Farm Animal, Wildlife and Food Safety Law . 2017.