Vol. 5, No.2, 76-85,  Januari  2009

## Analisis Kestabilan Model Populasi Mangsa Pemangsa dengan Pemanenan dan Waktu Tunda Melalui Pendekatan Bentuk
dt dx t x t x      ) ( ) (

## Syamsuddin Toaha †

## Abstrak
Dalam  tulisan  ini  dibahas  model  populasi  mangsa  pemangsa  dengan  pemanenan  usaha konstan  dan  waktu  tunda.  Bentuk ) (   t x yang  merupakan  pengaruh  waktu  tunda  pada model  didekati  oleh  bentuk dt dx t x   ) ( dengan  menggunakan  deret  Taylor.  Dengan pendekatan  itu,  model  tidak  lagi  mengandung  bentuk ) (   t x tetapi  tetap  bergantung  pada waktu  tunda.    Model  yang  dihasilkan  selanjutnya  dianalisis  kestabilannya  dengan menggunakan  nilai  eigen.    Untuk  model  tanpa  pemanenan,  diperoleh  suatu  titik keseimbangan yang positif dan stabil untuk
0   . Untuk 0   , titik keseimbangan mungkin tetap  stabil  dan  mungkin  juga  menjadi  tidak  stabil.    Untuk  model  yang  melibatkan pemanenan,  mungkin  diperoleh  suatu  titik  keseimbangan  yang  positif  asalkan  nilai  usaha pemanenan  dikontrol  dan  titik  keseimbangan  yang  terhasil  juga  stabil  untuk 0   . Sementara  jika 0   ,  titik  keseimbangan  mungkin  tetap  stabil  dan  mungkin  juga  menjadi tidak stabil. Keberadaan waktu tunda pada model mangsa pemangsa mungkin mempengaruhi kestabilan titik keseimbangan yang stabil menjadi tidak stabil.
Kata Kunci : Model mangsa pemangsa, pemanenan, waktu tunda,  nilai eigen.

## 1. Pendahuluan
Model  populasi  mangsa–pemangsa  adalah  salah  satu  model  yang  telah  dianalisis  secara meluaskan oleh banyak peneliti. Ini disebabkan karena banyak fenomena dalam interaksi populasi yang bersifat predasi (mangsa–pemangsa). Dalam penelitian ini model laju pertumbuhan populasi mangsa–pemangsa  yang  merujuk  kepada  model  Lotka–Volterra  dikembangkan  dengan mempertimbangkan  waktu  tunda  dan  pemanenan.  Waktu  tunda  dilibatkan  dalam  laju pertumbuhan  populasi  karena  pada  kenyataannya  laju  pertumbuhan  populasi  tidak  hanya bergantung  pada  jumlah  populasi  pada  saat  sekarang  tetapi  juga  bergantung  kepada  jumlah populasi sebelumnya. Waktu tunda dimasukkan ke dalam model untuk membuat model populasi menjadi lebih akurat yang lebih mendekati fenomena sebenarnya.
Dinamika  populasi  disamping  memberikan  konsekuensi  biologi  dan  ekologi,  juga memberikan  konsekuensi  ekonomi.  Populasi  ikan,  sebagai  contoh,  diharapkan  secara  ekologi tidak akan punah dan dari aspek ekonomi populasi diharapkan dapat memberikan banyak manfaat kepada  manusia.  Dengan  demikian  model  populasi  yang  mana  populasinya  bernilai  ekonomi seharusnya melibatkan fungsi pemanenan dalam model pertumbuhannya.
Dengan berbagai alasan, waktu tunda telah disetujui dan diterima untuk dipertimbangkan dalam berbagai model laju pertumbuhan populasi, lihat Bellman dan Cooke (1963), Hale (1977),
† Jurusan Matematika FMIPA Universitas Hasanuddin Makassar

## Syamsuddin Toaha
Cushing (1977), Gopalsamy (1992), Kuang (1993) dan MacDonald (1978) untuk sistem biologi dengan waktu tunda secara umum.
Persamaan  differensial  tunda  secara  umum  memperlihatkan  dinamika  yang  lebih kompleks  dibandingkan  dengan  persamaan  differensial  biasa  karena  waktu  tunda  dapat menyebabkan suatu titik keseimbangan yang stabil menjadi tidak stabil dan menyebabkan jumlah populasi  ikut  terpengaruh.  Waktu  tunda  yang  disebabkan  oleh  kehamilan  adalah  suatu  contoh yang  sederhana.  Secara  umum  konsumsi  mangsa  oleh  pemangsa  selama  waktu  yang  lalu menentukan  laju  perubahan  pemangsa  pada  waktu  sekarang.  Dengan  demikian,  model  interaksi suatu  populasi  dengan  populasi  lainnya  yang  lebih  realistik  seharusnya  mempertimbangkan pengaruh waktu tunda, Ho dan Ou (2002).
Dalam  tulisan  ini  yang  membahas  model  populasi  mangsa–pemangsa  yang  melibatkan waktu  tunda  dan  pemanenan,  akan  dianalisis syarat  kewujudan  titik  keseimbangan  yang  positif. Selanjutnya  digunakan  suatu  pengembangan  metode  dalam  menganalisis  pengaruh  waktu  tunda terhadap  kestabilan  titik  keseimbangan  model.  Tulisan  ini  bertujuan  untuk  mengontrol  nilai parameter, waktu tunda dan taraf pemanenan sehingga populasi mangsa dan pemangsa tetap akan wujud meskipun kedua populasi dieksploitasi dan dikenakan waktu tunda.
Metode  yang  digunakan  dalam  penelitian  ini  adalah  menghapuskan  bentuk  jumlah populasi  yang  bergantung  kepada  waktu  tunda  dengan  menggunakan  suatu  pendekatan  melalui deret  Taylor.  Selanjutnya  dilakukan  pelinearan  di  sekitar  titik  keseimbangan  dengan menggunakan  matriks  Jacobian.  Untuk  penentuan  kestabilan  titik  keseimbangannya  digunakan metode nilai eigen.

## 2. Review Model  Populasi Mangsa-Pemangsa dan Pemanenan
Model  populasi  mangsa  pemangsa  yang  berdasarkan  kepada  model  Lotka–Volterra adalah  salah  satu  model  yang  sangat  populer  dalam  matematika  biologi.  Luckinbill  (1973) menunjukkan  bahwa  populasi  mangsa  dan  pemangsa  dapat  hidup  bersama  dengan  cara mengurangi  frekuensi  interaksinya.  Danca et  al. (1997)  menganalisis  suatu  model  mangsa- pemangsa  dengan  menggunakan  metode  analitik  dan  numerik.  Mereka  dapat  menentukan  range nilai parameter dimana sistem mempunyai keadaan stasioner atau perilaku yang kacau.
Brauer  dan  Soudack  (1979a,  1979b  dan  1981)  mempertimbangkan  beberapa  model mangsa–pemangsa  yang  melibatkan  masalah  pemanenan  dengan  laju  konstan.  Mereka menganalisis  prilaku  global  model  mangsa–pemangsa  dan  membuat  klasifikasi  beberapa kemungkinan  dan  juga  menentukan  domain  atraksi  dari  trayektori.  Mereka  juga  menemukan bahwa  pada  kondisi  tertentu  model  dengan  pemanenan  adalah  stabil.  Ketika  model  mangsa- pemangsa dalam kondisi stabil asimptotik dan kestabilannya lemah maka limit cycle asimptotik kemungkinan akan wujud, Jeffries (1974).
Kar dan Chauduri (2004) mengkaji model mangsa–pemangsa yang juga berdasarkan pada model  Lotka–Volterra  dengan  pemanenan.  Mereka  menelaah  kemungkinan  kewujudan ekuilibrium  bionomik  dan  pemanenan  optimal.  Pengaruh  pemanenan  dengan  kuota  konstan  dan usaha konstan telah dikaji oleh Holmberg (1995) dan hasil yang diperoleh menunjukkan bahwa tangkapan  dengan  kouta  konstan  dapat  mengarah  kepada  keadaan  osilasi  dan  kacau  dan  dapat juga  meningkatkan  resiko  eksploitasi  yang  berlebihan.  Kestabilan  pemanenan  pada  model  yang melibatkan  dua  mangsa–satu  pemangsa  yang  berdasarkan  kepada  model  Lotka–Volterra  dalam Azar (1995) dimana populasi pemangsa dipanen juga telah dianalisis. Hasil analisis menunjukkan bahwa penggunaan fungsi pemanenan konstan membuat model lebih stabil.
Suatu  model  mangsa–pemangsa  dengan  tipe  Holling  yang  menggunakan  usaha pemanenan  sebagai  suatu  kontrol  disajikan  oleh  Srinivasu et  al .  (2001).  Mereka  menunjukkan bahwa  dengan  pemanenan,  ada  kemungkinan  untuk  menghentikan  perilaku  siklik  sistem  dan mengantarkan  suatu  limit  cycle  yang  stabil  global  dalam  sistem.  Hogarth et  al. (1992)

## Syamsuddin Toaha
mempertimbangkan  pemanenan  dengan  hasil  konstan  pada  sistem  satu  mangsa–satu  pemangsa dimana  mangsa  dan  pemangsa  dipanen  dengan  hasil  konstan  dan  hasil  analisis  menunjukkan bahwa kestabilan pada hasil maksimum yang berkelanjutan itu dibuktikan.
Matsuda  dan  Abrams  (2004)  menyelidiki  pengaruh  pada  ukuran  populasi  dan  hasil pemanenan pada level yang berbeda pada suatu pemangsa dalam sistem mangsa–pemangsa. Hasil analisis menunjukkan bahwa jumlah pemangsa mungkin meningkat dengan meningkatnya usaha penangkapan, meskipun ketika mangsa digolongkan memberikan suatu pengaruh positif pada laju pertumbuhannya ataupun ketika mangsa dieksploitasi berlebih oleh pemangsa.
Song  dan  Chen  (2002)  mempertimbangkan  esploitasi  pada  populasi  mangsa–pemangsa dengan  tingkat  struktur  dan  pemanenan  untuk  mangsa  dan  hasil  menunjukkan  bahwa  titik keseimbangan nonnegatif stabil global dan asimptotik  pada suatu kondisi tertentu. Di samping itu wujud  suatu  ambang  pemanenan  pada  populasi  mangsa  dan  waktu  tunda  optimal memaksimumkan jumlah total populasi.
Sistem  mangsa-pemangsa  telah  banyak  dikaji  penulis,  misalnya  oleh  Rinaldi  dan Muratori  (1992),  Kuang  dan  Beretta  (1998)  dan  Fan  dan  Kuang  (2004).  Sementara  model mangsa–pemangsa dengan pemanenan telah dikaji oleh Kitabatake (1987), Jonzen et al . (2001), Toaha  dan  Budin  (2005)  dan  Toaha et  al .  2007,  serta  model  mangsa–pemangsa  dengan  waktu tunda telah ditelaah Xu et al . (2004).

## 3. Model Mangsa Pemangsa dengan Pemanenan dan Waktu Tunda
Dalam penelitian ini kita akan mempertimbangkan suatu model mangsa–pemangsa yang berdasarkan  kepada  model  Lotka–Volterra  dengan  satu  populasi  mangsa  dan  satu  popoluasi pemangsa. Misalkan ) ( t x x  dan ) ( t y y  menyatakan ukuran populasi mangsa dan pemangsa pada  waktu t ,  maka  model  laju  pertumbuhan  populasi  mangsa  dan  pemangsa  diberikan  oleh sistem persamaan differensial autonomi
. 1 xy cy dt dy xy K x rx dt dx              
(1)
Model (1) meliputi parameter K , yaitu kapasitas bawaan populasi mangsa pada kondisi tidak ada interaksi  dengan  populasi  pemangsa.  Parameter r adalah  laju  pertumbuhan  intrinsik  populasi mangsa, c adalah  laju  kematian  populasi  pemangsa,  mengukur  laju  konsumsi  mangsa  oleh pemangsa,  mengukur konversi mangsa yang dikonsumsi ke dalam laju reproduksi pemangsa. Bermula  dari  model  logistik  waktu  tunda  oleh  Hutchinson,  May  (1974)  mengusulkan model mangsa–pemangsa dengan waktu tunda
         
        . 1 t y t x t cy dt t dy t y t x K t x t rx dt t dx                 (2)
Model  (2)  meliputi  satu  waktu  tunda  dan  laju  pertumbuhan  populasi  mangsa  tidak  hanya bergantung kepada ukuran populasi sekarang tetapi juga bergantung kepada ukuran populasi pada

## Syamsuddin Toaha
beberapa waktu sebelumnya. Jika periode persiapan populasi mangsa untuk bereproduksi adalah  , maka laju pertumbuhan per kapita harus memuat suatu waktu tunda  . Model (2) merupakan pengembangan dari model (1) dimana laju pertumbuhan populasi mangsa tidak hanya bergantung pada  ukuran  populasi  sekarang  tetapi  juga  bergantung  pada  ukuran  populasi  yang  waktu sebelumnya.
Model  (1)  kemudian  dikembangkan  dengan  pertimbangan  bahwa  kedua  populasi mempunyai nilai jual, sebagai contoh populasi ikan. Disini, kedua populasi dipanen dengan usaha konstan. Pengaruh pemanenan ini tidak serta merta mempengaruhi laju perubahan populasi, tetapi perlu  beberapa  waktu  untuk  memberikan  efek  kepada  laju  pertumbuhan  populasi.  Dengan demikian laju tangkapan proporsional dengan ukuran populasi pada beberapa waktu sebelumnya. Model  populasi  mangsa–pemangsa  yang  terhasil  dengan  waktu  tunda  dalam  bentuk  pemanenan adalah
                  . ) ( ) ( 1                     t y E q t y t x t cy dt t dy t x E q t y t x K t x t rx dt t dx y y x x (3)
Konstanta q x dan q y adalah koefisien ketertangkapan populasi mangsa dan pemangsa manakala E x dan E y adalah usaha pemanenan untuk populasi mangsa dan pemangsa. Penglibatan waktu tunda dalam  bentuk  laju  pemanenan  juga  telah  dipertimbangkan  oleh  Kar  (2003).  Untuk  keperluan analisis model diasumsikan 1   y x q q .
Dalam  kasus  dimana  populasi  mangsa  dan  pemangsa  ditangkap  dengan  menggunakan alat  penangkapan  yang  sama,  maka  sangat  beralasan  untuk  mengasumsikan  bahwa  usaha penangkapan yang digunakan adalah sama, yaitu y x E E E   .
Dari pengembangan model populasi mangsa–pemangsa, ada beberapa permasalahan yang akan  dianalisis  dan  diselesaikan,  yaitu  (1)  menentukan  syarat  kewujudan  titik  keseimbangan model  yang  positif,  (2)  menentukan  nilai  waktu  tunda  sehingga  titik  keseimbangan  yang  stabil tetap stabil, (3) menentukan interval untuk nilai usaha pemanenan sehingga tetap diperoleh titik keseimbangan yang positif dan stabil, (4) mengontrol nilai waktu tunda, parameter model, jumlah awal  populasi  dan jumlah tangkapan  sehingga  wujud  titik  keseimbangan  yang  positif  dan  stabil sehingga populasi tidak akan punah untuk selamanya.

## 4. Metodologi dan Analisis Kestabilan
Model  populasi  mangsa–pemangsa  akan  dianalisis  untuk  menentukan  kestabilan  titik keseimbangan  dan  menentukan  ambang  batas  waktu  tunda  sehingga  titik  keseimbangan  yang stabil tetap stabil. Sebelum penentuan kestabilan titik keseimbangan model, kita sediakan syarat perlu atau cukup untuk kewujudan titik keseimbangan positif karena model yang dianalisis adalah model populasi.
Metode  yang  digunakan  dalam  menganalisis  kestabilan  titik  keseimbangan  dan penentuan  ambang  batas  waktu  tunda  adalah  metode  pelinearan,  penghapusan  bentuk ) (   t x dengan  melakukan  pendekatan  dan  transformasi.  Dari  model  linear  yang  terhasil,  penentuan kestabilan titik keseimbangan dan ambang batas waktu tunda dapat dilakukan.
Langkah-langkah dalam metode ini secara berurut adalah sebagai berikut :
1. Penentuan titik  keseimbangan  model  yang  positif  untuk  model  nonlinear  yang  mengandung waktu tunda.

## Syamsuddin Toaha
2. Menentukan  suatu  syarat  sehingga  titik  keseimbangan  yang  positif  itu  stabil  jika  tidak  ada waktu tunda.
3. Menghapuskan bentuk ) (   t x dengan pendekatan dt dx t x t x      ) ( ) ( . 4. Menggantikan ) (   t x dengan
dt dx t x   ) ( pada  model  sehingga  diperoleh  model  tanpa bentuk ) (   t x tetapi masih melibatkan waktu tunda  .
5. Melinearkan  model  di  sekitar  titik  keseimbangan  positif  dengan  menggunakan  matriks Jacobian sehingga terhasil model linear dengan titik keseimbangan ) 0 , 0 ( .
6. Dari  model  linear  selanjutnya  ditentukan  syarat  kestabilan  titik  keseimbangan  melalui persamaan karakteristik dan nilai eigen dari matriks Jacobian.
7. Ambang batas nilai waktu tunda ditentukan dengan memberikan suatu syarat pada persamaan karakteristik sehingga semua nilai eigen mempunyai bagian riil yang negatif.
Metode  di  atas  dapat  diaplikasikan  untuk  menganalisis  kestabilan  titik  keseimbangan model  nonlinear  jika  titik  keseimbangan  itu  stabil  untuk 0   .  Selanjutnya,  karena  nilai  eigen yang terhasil merupakan suatu fungsi yang kontinu terhadap  maka nilai eigen yang bergantung kepada  juga masih mempunyai bagian riil yang positif untuk suatu waktu  yang relatif kecil. Dengan  menggunakan  pendekatan  yang  digunakan  dalam  metode  di  atas,  penentuan  ambang batas  nilai  waktu  tunda  dapat  dilakukan.  Jika  waktu  tunda  yang  diaplikasikan  dalam  model melebihi nilai ambang batas waktu tunda tersebut, maka titik keseimbangan menjadi tidak stabil yaitu pada saat terdapat nilai eigen dengan bagian riil positif, lihat Willems (1970).
Model (1) mempunyai titik keseimbangan positif        K c K r c    ) ( , jika 0    c K dan
titik  keseimbangan  ini  stabil  asimptotik  global,  Ho  dan  Ou  (2002).  Dengan  mengenakan pendekatan dt dx t x t x      ) ( ) ( pada bentuk waktu tunda untuk model (2) diperoleh
   
            , ) ( ) ( ) ( t y t x t cy dt t dy t y t x dt t dx t x t x K r t rx dt t dx                 yang dapat ditulis sebagai                 . ) ( 1 ) ( 2 t y t x t cy dt t dy t x K r t y t x t x K r t rx dt t dx
                (4)
Titik  keseimbangan  positif  untuk  model  (4)  sama  dengan  titik  keseimbangan  model  (1) ataupun  model  (2)  yaitu        K c K r c    ) ( , ,  namun  model  (4)  mensyaratkan  bahwa  r K t x  ) ( .
Matriks Jacobian untuk model (4) adalah
                    x c y K x r Kx K x r yK x r xrK rK J        ) ( ) ( 2 2 2 2 2 2 .
Matriks  Jacobian  yang  bersesuaian  dengan  titik  keseimbangan        K c K r c    ) ( , adalah
       0 1 R Q P J dengan K rc rc P     ,
K rc Kc Q      dan K c K r R   ) (   .  Nilai  eigen  dari
matriks  Jacobian 1 J adalah 2 4 2 2 , 1 QR P P     .  Untuk  analisis  kestabilan  titik keseimbangan ini, akan ditinjau 2 kasus.
Kasus 1:   r K c  atau rc K    .
Untuk rc K    maka diperoleh P dan Q keduanya bernilai negatif yang mengakibatkan bagian  real  dari  kedua  nilai  eigen  bernilai  negatif.  Dengan  demikian  titik  keseimbangan
       K c K r c    ) ( , stabil asimptotik jika rc K     0 . Mudah  dicek  bahwa 0 4 2   QR P jika  dan  hanya  jika     dimana
) ( 4 ) ( 4 c K rc rc c K K          .  Dengan  demikian  titik  keseimbangan   
    K c K r c    ) ( , stabil
asimptotik  dan  merupakan  node  jika   rc K   
   , 0 max dan  stabil  spiral  jika      0
sekiranya   bernilai positif.
Kasus 2:   r K c  atau rc K    .
Untuk rc K    maka diperoleh P dan Q keduanya bernilai positif yang mengakibatkan kedua  nilai  eigen  bernilai  real  dan  berbeda  tanda.  Dengan  demikian  titik  keseimbangan
       K c K r c    ) ( , tidak stabil yang merupakan titik pelana.
Dengan mengasumsikan 1   y x q q dan y x E E
E   , maka model (3) dapat ditulis sebagai                 . ) ( ) ( ) ( 2 t y t x t Ey t cy dt t dy t y t x t x K r t Ex t rx dt t dx
   
          (5)

## Syamsuddin Toaha
Dengan  mengenakan  pendekatan
dt t dx t x t x ) ( ) ( ) (      pada  bentuk  tundaan  waktu  untuk model (5) diperoleh
       
        , ) ( ) ( ) ( ) ( ) ( 2 t y t x dt t dy t y E t cy dt t dy t y t x t x K r dt t dx t x E t rx dt t dx                           yang dapat ditulis sebagai                     , 1 1 ) (
1 2 1     E t y t x t y c dt t dy E t y t x t x K r t x r dt t dx
        (6)
dengan E r r   1 dan E c c   1 serta diasumsikan bahwa 0   E r .
Dengan  menganggap  bahwa E 1   ,  model  (6)  mempunyai  titik  keseimbangan  positif
       K rc K r c    ) , 1 1 1 apabila 0 1 1   rc K r  dan  hal  ini  terjadi  jika 1 E E  dengan
dan  diasumsikan  pula  bahwa 0   c K  .  Matriks  Jacobian  untuk  model  (6) adalah
                      ) 1 ( ) 1 ( ) 1 ( ) 1 ( 2 2 1 1         E c x E y E x E K yK rx K r J .
Matriks Jacobian yang bersesuaian dengan titik keseimbangan        K rc K r c    ) , 1 1 1 adalah        0 3 1 1 1 R Q P J
dengan ) 1 ( 1 1     E K rc P , ) 1 ( 1 1      E K c Q dan ) 1 ( ) ( 1 1 1        E K rc K r R .  Nilai  eigen  dari
matriks  Jacobian 3 J adalah 2 4 1 1 2 1 1 2 , 1 R Q P P     . Untuk  analisis  kestabilan  titik keseimbangan ini, akan ditinjau 2 kasus.
) ( ) ( 1 r K c K r E     
Kasus 1 . 0 1    E atau
E 1   . Untuk E 1   diketahui bahwa 1 P dan 1 Q bernilai negatif sementara 1 R bernilai positif. Dengan  demikian  diperoleh  bagian  real  dari  nilai  eigen  matriks 3 J keduanya  bernilai  negatif yang berarti bahwa titik keseimbangan        K rc K r c    ) , 1 1 1 stabil asimptotik. Mudah  dicek  bahwa 0 4 1 1 2 1   R Q P jika  dan  hanya  jika
2 E E  dengan ) ( 4 ) ( 4 2 c K rc rc c K K E        .  Dengan  demikian  titik  keseimbangan        K rc K r c    ) , 1 1 1 stabil
asimptotik  dan  merupakan node jika  memenuhi E 1   dan 2 0 E E   .  Dapat  dengan  mudah dibuktikan bahwa 1 2 E E  . Sementara jika E 1   dan 1 2 E E E   maka titik keseimbangan
       K rc K r c    ) , 1 1 1 stabil spiral.
Kasus 2 . 0 1    E atau E 1   . Untuk E 1   diketahui  bahwa 1 P dan 1 Q keduanya  bernilai  positif  dan 1 R juga bernilai positif. Dengan demikian diperoleh nilai eigen matriks 3 J keduanya bernilai real dengan tanda  yang  berbeda,  positif  dan  negatif,  yang  berarti  bahwa  titik  keseimbangan
       K rc K r c    ) , 1 1 1 tidak stabil, ia adalah titik pelana.

## Daftar Pustaka
1. Azar,  C.  1995.  Long-term  environmental  problems  economic  measures  and  physical indicators.  http://www2.lib.chalmers.se/cth/diss/doc/9596/Azar  Christian.html [  20 Agustus 2004].
2. Bellman, R. dan Cooke, K.L. 1963. Differential difference equations . New York: Academic Press.
3. Brauer,  F.  dan  Soudack,  A.C.  1979a.  Stability  regions  and  transition  phenomena  for harvested predator-prey systems. Journal Math. Biology 7:319-337.
4. Brauer,  F.  dan  Soudack,  A.C.  1979b.  Stability  regions  in  predator-prey  systems  with constant-rate prey harvesting. Journal Math. Biology 8:55-71.
5. Brauer, F. dan Soudack, A.C. 1981. Coexistence properties of some predator-prey systems under constant rate harvesting and stocking. Journal Math. Biology 12:101-114.

## Syamsuddin Toaha
6. Cushing,  J.M.  1977. Integrodifferential  equations  and  delay  models  in  population dynamics . Heidelberg: Springer-Verlag.
7. Danca, M., Codreanu, S. dan Bako, B. 1997. Detailed analysis of a nonlinear prey-predator model. Journal of Biological Physics 23:11-20.
8. Fan, M. dan Kuang, Y. 2004. Dynamics of a nonautonomous predator-prey system with the Beddington-DeAngelis functional response. J. Math. Anal. Appl. 295:15-39.
9. Gopalsamy, K. 1992. Stability and oscillations in delay differential equations of population dynamics . Netherlands: Kluwer Academic Publishers.
10. Hale, J.K. 1977. Theory of functional differential equations . Heidelberg: Springer-Verlag.
11. Ho, C.P. dan Ou, Y.L. 2002. Influence of time delay on local stability for a predator-prey system. Journal of Tunghai Science 4:47-62.
12. Hogarth, W.L., Norbury, J., Cunning, I. dan Sommers, K. 1992. Stability of a predator-prey model with harvesting. Ecological Modelling 62:83-106.
13. Holmberg,  J.  1995.  Socio-ecological  principles  and  indicators  for  sustainability,  PhD Thesis, Goteborg University, Sweden.
14. Jeffries, C. 1974. Probabilistic limit cycles. In Lecture notes in biomathematics (2nd ed.), ed. S. Levin, pp 123-131. New York: Spinger-Verlag.
15. Jonzen,  N.,  Lundberg,  P.  dan  Gardmark,  A.  2001.  Harvesting  spatially  distributed populations. J. Wildl. Biol . 7:197-203.
16. Kar,  T.K.  2003.  Selective  harvesting  in  a  prey-predator  fishery  with  time  delay. Mathematical and Computer Modelling 38:449-458.
17. Kar, T.K. dan Chaudhuri, K.S. 2004. Harvesting in a two-prey one predator fishery: A bioeconomic model. J. ANZIAM 45:443-456.
18. Kitabatake, Y. 1987. A  dynamic predator-prey  model for the utilization of fishery resources:  A case of trawling in  lake  Kasumigaura.  In Mathematical modelling  of environmental  and  ecological  systems ,  ed.  J.B.  Shukla,  T.G.  Hallam,  and  V. Capasso, pp 233-253. Elsevier Science Publishers.
19. Kuang,  Y.  1993. Delay  differential  equations  with  application  in  population dynamics . New York: Academic Press.
20. Kuang,  Y.  dan  Beretta,  E.  1998.  Global  qualitative  analysis  of  a  ratio-dependent predator-prey system. J. Math. Biol . 36:389-406.
21. Luckinbill,  L.  S.  1973.  Coexistence  in  laboratory  populations  of  paramecium aurelia and its predator didinium nasutum. Journal of Ecology 54(6):1320-1327.

## Syamsuddin Toaha
22. MacDonald, N. 1978. Time lag in biological models . Heidelberg: Springer-Verlag.
23. Matsuda,  H.  dan  Abrams,  P.A.  2004.  Effects  of  predators-prey  interaction  and adaptive change on sustainable yield. Can. J. Fish. Aquat. Sci./J. Can. Sci. Halieut. Aquat. 61(2):175-184.
24. May,  R.M.  1974. Stability  and c omplexity  of  model  ecosystems .  Princeton,  New Jersey: Princeton University Press.
25. Rinaldi,  S.  dan  Muratori,  S.  1992.  Slow-fast  limit  cycles  in  predator-prey  models. Ecological Modelling 61:287-308.
26. Song,  X.Y.  dan  Chen,  L.S.  2002.  Optimal  harvesting  and  stability  for  a  predator- prey system with stage structure. Acta Mathematicae Applicatae Sinica 18(3):423- 430.
27. Srinivasu,  P.D.,  Ismail,  S.  dan  Naidu,  C.R.  2001.  Global  dynamics  and controllability of a harvested prey-predator system. J. Biological Systems 9(1):67- 79.
28. Toaha,  S.  dan  Budin,  H.  2005.  Stability  analysis  of  one  prey  –  two  predators population model under constant quota of harvesting . In Proceedings of The second International Conference on Research and Education in Mathematics (ICREM 2) , ed. U.D. Bekbaev, pp 365-374. Serdang, Malaysia.
29. Toaha,  S.,  Hassan,  M.A.,  Ismail,  F.  dan  June,  L.W.  2007.  Stability  analysis  and maximum  profit  of  one  prey-  two  predators  model  under  constant  effort  of harvesting, Malaysian Journal of Science 26:43-51.
30. Willems, J.L. 1970. Stability theory of dynamical systems . London: Thomas Nelson & Sons.
31. Xu,  R.,  Chaplain,  M.A.J.  dan  Davidson,  F.A.  2004.  Periodic  solution  for  a  three- species  Lotka-Voterra  food-chain  model  with  time  delays. Mathematical  and Computer Modelling 40:823-837.