Bifurkasi Hopf pada Model Lotka-Volterra Orde-Fraksional dengan Efek Allee Aditif pada Predator
Hasan S. Panigoro 1,3* , Dian Savitri 2,3
1 Jurusan Matematika, Universitas Negeri Gorontalo, Gorontalo 96119, Indonesia
2 Jurusan Matematika, Universitas Negeri Surabaya, Surabaya 60231, Indonesia
3 Jurusan Matematika, Universitas Brawijaya, Malang 65145, Indonesia
* Penulis Korespondensi: hspanigoro@ung.ac.id

## Jambura Journal of Biomathematics
Jambura J. Biomath. Volume 1, Issue 1, pp. 16-24, June 2020
Journal Homepage: http://ejurnal.ung.ac.id/index.php/JJBM/index https://doi.org/10.34312/jjbm.v1i1.6908
Received: 12 June 2020,
Accepted: 17 June 2020, Published Online: 29 June 2020

## Abstrak
Artikel ini bertujuan untuk mempelajari dinamika dari model predator-prey Lotka-Volterra dengan adanya efek Allee pada predator . Berdasarkan kondisi biologis, dipilih turunan Caputo orde-fraksional sebagai operatornya. Analisis terhadap model diawali dengan mengidentifikasi eksistensi, ketunggalan, dan kepositifan solusi. Selanjutnya, dilakukan investigasi terhadap eksistensi titik kesetimbangan dan kestabilannya. Diperoleh titik kepunahan kedua populasi yang selalu tidak stabil tipe pelana, dan titik eksistensi kedua populasi yang secara kondisional stabil baik lokal maupun global. Salah satu fenomena yang menarik pada model ini adalah perubahan kestabilan pada titik eksistensi kedua populasi yang diakibatkan oleh perubahan dari orde turunannya yang disebut dengan bifurkasi Hopf. Pada akhirnya, diberikan simulasi numerik untuk melakukan validasi terhadap hasil teoritis.
Kata Kunci: Lotka-Volterra; Efek Allee Aditif; Bifurkasi Hopf; Orde-Fraksional

## Abstract
This article aims to study the dynamics of a Lotka-Volterra predator-prey model with Allee effect in predator. According to the biological condition, the Caputo fractional-order derivative is chosen as its operator. The analysis is started by identifying the existence, uniqueness, and non-negativity of the solution. Furthermore, the existence of equilibrium points and their stability is investigated. It has shown that the model has two equilibrium points namely both populations extinction point which is always a saddle point, and a conditionally stable co-existence point, both locally and globally. One of the interesting phenomena is the occurrence of Hopf bifurcation driven by the order of derivative. Finally, the numerical simulations are given to validate previous theoretical results.
Keywords: Lotka-Volterra; Additive Allee Effect; Hopf Bifurcation; Fractional-Order

## 1. Pendahuluan
Fenomena efek Allee merujuk pada penurunan laju pertumbuhan suatu populasi pada saat kepadatan populasinya rendah [1–4] . Kondisi ini pertama kali dipelajari oleh Warder Cyle Allee pada tahun 1931 [5] . Penurunan laju pertumbuhan ini disebabkan oleh beberapa kondisi biologis seperti kompetisi intraspesifik, kesulitan mendapatkan pasangan, dan turunnya kewaspadaan terhadap predator alaminya, lihat ref. [6–9] . Efek Allee dapat terjadi pula pada interaksi populasi, baik pada prey maupun predator , sehingga mengakibatkan perubahan perilaku populasi di bidang ekologi. Stephens dkk [10] melakukan konstruksi model matematika efek Allee aditif serta konsekuensi dalam kelangsungan hidup suatu populasi khususnya terhadap konservasi. Beberapa penelitian menunjukkan bahwa efek Allee dapat ditemukan dalam suatu sistem ekologi seperti pada model predator-prey [11–14] . Dengan demikian adanya efek Allee memberikan perubahan terhadap dinamika dari model apabila konstanta Alleenya divariasikan. Salah satu model predator-prey yang cukup dikenal adalah model klasik Lotka-Volterra [15] . Model ini didefinisikan oleh
dx dt = ax − bxy , dy dt = cxy − dy , (1)
dengan x ( t ) adalah kepadatan populasi prey dan y ( t ) adalah kepadatan populasi predator . Parameter a , b , c , dan d masing-masing merupakan laju pertumbuhan prey , laju pemangsaan, laju konversi pemangsaan terhadap pertumbuhan predator , dan laju kematian alami predator . Secara biologis, efek Allee dapat terjadi pada
prey maupun predator . Pada artikel ini, diasumsikan bahwa populasi yang mengalami efek Allee adalah predator . Secara matematis, efek Allee yang digunakan adalah efek Allee aditif seperti pada ref. [1] dan [6] . Dengan demikian diperoleh model sebagai berikut.
dx dt = ax − bxy , dy dt = cxy − dy − my y + n , (2)
dimana m dan n adalah konstanta yang menyatakan efek Allee. Beberapa penelitian menunjukkan bahwa untuk mendapatkan model yang lebih realistis, digunakan model dengan turunan orde-fraksional seperti pada ref. [16–19] . Didefinisikan turunan orde-fraksional Caputo sebagai berikut.
Definisi 1. [20] Misalkan 0 < α ≤ 1. Turunan fraksional Caputo orde- α didefinisikan oleh
C D α t u ( t ) = 1 Γ ( 1 − α ) Z t 0 ( t − s ) − α u 0 ( τ ) d τ , (3)
dimana t ≥ 0, u ∈ C n ([ 0, + ∞ ) , R ) , dan Γ ( · ) adalah fungsi Gamma.
Dengan mengganti ruas kiri dari model (2) dengan turunan fraksional Caputo seperti pada Definisi 1 , diperoleh model orde-fraksional sebagai berikut.
C D α t x = ax − bxy = F 1
C D α t y = cxy − dy − my y + n = F 2 (4)
Model (4) merupakan objek dan kajian utama pada artikel ini. Untuk mendukung tahapan analisis, terlebih dahulu diberikan beberapa teorema dan lema. Setelah itu, diberikan hasil-hasil analisis berupa eksistensi, ketunggalan, dan kepositifan solusi, kestabilan lokal dan kestabilan global dan bifurkasi Hopf. Di akhir pembahasan, diberikan simulasi numerik dari model (4) yang bersesuaian dengan hasil analitiknya.

## 2. Metode
Dalam mempelajari dinamika dari model (4) , ada beberapa tahapan yang dilakukan yaitu:
1. Mengidentifikasi eksistensi, ketunggalan, dan kepositifan solusi model (4)
2. Mengidentifikasi eksistensi dari titik kesetimbangan, kestabilan lokal, dan kestabilan global dari model (4) . 3. Mempelajari terjadinya bifurkasi Hopf pada model (4) .
4. Memberikan simulasi numerik dan interpretasi biologis atas kajian teoritis sebelumnya.
Untuk dapat menyelesaikan tahapan-tahapan tersebut, diberikan beberapa teorema dan lemma yang dapat mendukung proses analisis yaitu sebagai berikut.
Lema 2. [21] Misalkan α ∈ ( 0, 1 ] , u ( t ) ∈ C [ p , q ] dan C D α t u ( t ) ∈ C [ p , q ] . Jika C D α t u ( t ) ≥ 0, ∀ t ∈ ( p , q ) , maka u ( t ) adalah fungsi tak-turun untuk setiap t ∈ [ p , q ] . Jika C D α t u ( t ) ≤ 0, ∀ t ∈ ( p , q ) , maka u ( t ) adalah fungsi tak-naik untuk setiap t ∈ [ p , q ] .
Teorema 3. (Kondisi Matignon [20, 22] ) Misalkan α ∈ ( 0, 1 ] . Pandang sistem persamaan diferensial fraksional orde- α
C D α t ~ x = ~ f ( ~ x ) , ~ x ( 0 ) = ~ x 0 . (5)
~ x ∗ adalah titik kesetimbangan dari pers. (5) jika memenuhi ~ f ( ~ x ∗ ) = 0 . Titik kesetimbangan ini stabil asimtotik lokal jika semua nilai eigen λ j dari matriks Jacobi J = ∂ f ∂ x di ~ x ∗ memenuhi | arg ( λ j ) | > απ 2 . Jika terdapat masing-masing paling sedikit satu nilai eigen yang memenuhi | arg ( λ k ) | > απ
2 dan | arg ( λ l ) | < απ 2 , k 6 = l, maka ~ x ∗ disebut titik pelana.
Lema 4. [23] Pandang sistem persamaan diferensial berikut
C D α t x ( t ) = g ( t , x ( t )) , t > 0, x ( 0 ) ≥ 0, α ∈ ( 0, 1 ] , (6)
dimana g : ( 0, ∞ ) × Φ → R n , Φ ⊆ R n . Sistem pers. (6) memiliki solusi tunggal di ( 0, ∞ ) × Φ jika g ( t , x ( t )) memenuhi kondisi Lipschitz lokal.
Lema 5. [24] Misalkan u ( t ) ∈ C ( R + ) , x ∗ ∈ R + , dengan turunan fraksional Caputo orde- α ada untuk α ∈ ( 0, 1 ] . Maka, untuk suatu t > 0
C D α t u ( t ) − u ∗ − u ∗ ln u ( t ) u ∗ ≤ 1 − u ∗ u ( t ) C D α t u ( t ) .
Lema 6. [25] ) Misalkan Ω merupakan himpunan tertutup dan terbatas dengan semua solusi dari sistem
C D α t x ( t ) = f ( x ( t )) , (7)
dengan nilai awal di Ω dan terus berada di Ω sepanjang waktu. Jika ∃ V ( x ) : Ω → R yang memiliki turunan parsial dan memenuhi C D α t V | ( 7 ) ≤ 0 . Misalkan E : = n x | C D α t V | ( 7 ) = 0 o dan M adalah himpunan invarian terbesar dari E, maka setiap solusi dari Ω akan menuju M saat t → ∞ .

## 3. Hasil dan Pembahasan
Sebelum mempelajari dinamika dari model (4) , akan ditunjukkan bahwa solusi dari model (4) ada dan tunggal, serta selalu bernilai positif untuk nilai awal tak negatif. Berdasarkan alasan biologis, maka seluruh solusi harusnya berada di R 2 + dengan
R 2 + : = { ( x , y ) : x ≥ 0, y ≥ 0, x , y ∈ R } (8)

## 3.1. Eksistensi, Ketunggalan dan Kepositifan Solusi
Pada tahapan ini akan ditunjukkan bahwa untuk setiap nilai awal di Ψ : = ( x , y ) ∈ R 2 : max {| x | , | y |} ≤ γ maka eksistensi dan ketunggalan solusi dari model (4) selalu terpenuhi.
Teorema 7. Untuk setiap nilai awal ( x 0 , y 0 ) ∈ R 2 + di Ψ , terdapat solusi tunggal di Ψ dari model (4) sepanjang t > 0 . bukti . Akan digunakan pendekatan yang sama seperti di ref. [26] . Pandang suatu pemetaan F ( X ) = ( F 1 ( X ) , F 2 ( X )) . Untuk suatu X = ( x , y ) , ¯ X = ( ¯ x , ¯ y ) , dengan X , ¯ X ∈ Ψ maka dari model (4) diperoleh k F ( X ) − F ( ¯ X ) k = | F 1 ( X ) − F 1 ( ¯ X ) | + | F 2 ( X ) − F 2 ( ¯ X ) | = | ax − bxy − ( a ¯ x − b ¯ x ¯ y ) | + cxy − dy − my y + n − c ¯ x ¯ y − d ¯ y − m ¯ y ¯ y + n = | a ( x − ¯ x ) − b ( xy − ¯ x ¯ y ) | + c ( xy − ¯ x ¯ y ) − d ( y − ¯ y ) − my y + n − m ¯ y ¯ y + n ≤ a | x − ¯ x | + ( b + c ) | xy − ¯ x ¯ y | + d | y − ¯ y | + m n | y − ¯ y | ≤ a | x − ¯ x | + ( b + c ) γ | y − ¯ y | + ( b + c ) γ | x − ¯ x | + d + m n | y − ¯ y | = ( a + ( b + c ) γ ) | x − ¯ x | + ( b + c ) γ + d + m n | y − ¯ y | ≤ L k X − ¯ X k
dengan L = max n a + ( b + c ) γ , ( b + c ) γ + d + m n o . Dengan demikian maka F ( X ) memenuhi kondisi Lipschitz. Berdasarkan Lema 4 maka untuk setiap nilai awal ( x 0 , y 0 ) ∈ R 2 + di Ψ , terdapat solusi tunggal di Ψ dari model (4) sepanjang t > 0 .
Teorema 8. Untuk semua nilai awal ( x 0 , y 0 ) ∈ R 2 + mengakibatkan semua solusi model (4) tak negatif.
bukti . Pertama-tama akan ditunjukkan bahwa dengan ( x 0 , y 0 ) ∈ R 2 + maka x ( t ) ≥ 0 sepanjang t → ∞ dengan menggunakan kontradiksi. Andaikan pernyataan tersebut tidak benar maka terdapat t 1 > 0 sehingga
   x ( t ) > 0, 0 ≤ t < t 1 x ( t 1 ) = 0, x ( t + 1 ) < 0. (9)
Dengan menggunakan persamaan (9) dan persamaan pertama di model (4) diperoleh
C D α t x ( t 1 ) | x ( t 1 ) = 0. (10)
Berdasarkan Lema 2 , diperoleh x ( t + 1 ) = 0, yang berarti kontradiksi dengan fakta bahwa x ( t + 1 ) < 0. Maka haruslah x ( t ) ≥ 0 sepanjang t → ∞ . Dengan mengunakan cara yang serupa maka dapat ditunjukkan bahwa y ( t ) ≥ 0 sepanjang t → ∞ .

## 3.2. Titik Kesetimbangan dan Kestabilannya
Titik kesetimbangan biologis diperoleh dengan mencari solusi positif dari sistem persamaan berikut.
[ a − by ] x = 0, cx − d − m y + n y = 0, (11)
sehingga diperoleh dua titik kesetimbangan yaitu E 0 = ( 0, 0 ) yang menunjukkan kepunahan kedua populasi dan E 1 = m a + bn + d b , a b yang menunjukkan eksistensi dari kedua populasi. Kedua titik ini merupakan titik kesetimbangan biologis yang berarti selalu terdefinisi di R 2 + . Selanjutnya akan dipelajari kestabilan lokal dan global dari kedua titik kesetimbangan tersebut.
Teorema 9. Titik E 0 = ( 0, 0 ) merupakan titik pelana.
bukti . Dengan mencari matriks Jacobi dari model (4) di titik E 0 = ( 0, 0 ) diperoleh
J ( E 0 ) = " a 0 0 − dn + m n # ,
yang memberikan nilai eigen λ 1 = a dan λ 2 = − dn + m n sehingga | arg ( λ 1 ) | = 0 < απ 2 dan | arg ( λ 2 ) | = π > απ 2 . Berdasarkan Teorema 3 , E 0 = ( 0, 0 ) merupakan titik pelana.
Berdasarkan Teorema 9 , populasi dari predator dan prey tidak akan pernah mengalami kepunahan. Eksistensi dari keduanya akan terjaga meskipun dengan kondisi awal yang cukup kecil.
Teorema 10. Misalkan
α ∗ = 2 π tan − 1   q 4 a ( bdn + ad + bm )( a + bn ) 3 − abm 2 abm  
(i) Untuk abm 2 ( a + bn ) 2 ! 2 ≥ a ( bdn + ad + bm ) a + bn , E 1 = m a + bn + d b , a b selalu tidak stabil.
(ii) Untuk abm 2 ( a + bn ) 2 ! 2 < a ( bdn + ad + bm ) a + bn , jika α ≤ α ∗ maka E 1 = m a + bn + d b , a b stabil asimtotik lokal dan jika α > α ∗ maka E 1 = m a + bn + d b , a b tidak stabil.
bukti . Matriks Jacobi dari model (4) di titik E 1 = m a + bn + d b , a b adalah
J ( E 0 ) =     0 − b ( bdn + ad + bm ) c ( a + bn ) ac b abm ( a + bn ) 2     ,
sehingga diperoleh persamaan karaketeristik λ 2 − abm
( a + bn ) 2 λ + a ( bdn + ad + bm ) a + bn = 0 yang memberikan nilai eigen
λ 1,2 = abm 2 ( a + bn ) 2 ± v u u t abm 2 ( a + bn ) 2 ! 2 − a ( bdn + ad + bm ) a + bn .
Jika abm 2 ( a + bn ) 2 ! 2 ≥ a ( bdn + ad + bm ) a + bn maka λ 1 λ 2 = a ( bdn + ad + bm ) a + bn > 0 dan λ 1 + λ 2 = abm ( a + bn ) 2 > 0 sehingga λ 1,2 > 0. Akibatnya | arg ( λ 1,2 ) | = 0 < απ 2 . Berdasarkan Teorema 3 , E 1 tidak stabil. Untuk abm
2 ( a + bn ) 2 ! 2 < a ( bdn + ad + bm ) a + bn , jika α ≤ α ∗ maka | arg ( λ 1,2 ) | > απ 2 . Berdasarkan Teorema 3 , E 1 stabil asimtotik lokal. Jika α > α ∗ maka | arg ( λ 1,2 ) | < απ 2 sehingga berdasarkan Teorema 3 , E 1 tidak stabil.
Perhatikan bahwa titik kesetimbangan E 1 selalu ada, sehingga berdasarkan Teorema 10 , diperoleh syarat biologis agar eksistensi dari kedua populasi terjaga dan konvergen ke titik kesetimbangan E 1 jika kondisi awalnya cukup dekat dengan E 1 . Syarat ini dianggap kurang kuat karena kondisi awal dari kedua populasi tidak selamanya cukup dekat dengan titik kesetimbangan E 1 . Oleh karena itu dubutuhkan syarat yang mengakibatkan kepadatan populasi selalu terjaga meskipun dengan kondisi awal seperti apapun. Hal ini ditunjukkan oleh teorema berikut.
Teorema 11. Titik E 1 = ( x ∗ , y ∗ ) = m a + bn + d b , a b stabil asimtotik global jika dn + m cn < x ∗ < d c .
bukti . Didefinisikan fungsi Lyapunov
V ( t ) = x − x ∗ − x ∗ ln x ( t ) x ∗ + b c y − y ∗ − y ∗ ln x ( t ) y ∗ (12)
Menurut Lema 5 diperoleh
C D α t V ( t ) ≤ ( x − x ∗ ) ( a − by ) + b c ( y − y ∗ ) cx − d − m y + n = ax − ax ∗ + bx ∗ y − by ∗ x − bdy c + bdy ∗ c − bmy c ( y + n ) + bmy ∗ c ( y + n ) = − ( by ∗ − a ) x − ax ∗ − bdy ∗ c − bd c − bx ∗ y − bmy c ( y + n ) + bmy ∗ c ( y + n ) ≤ − ( by ∗ − a ) x − ax ∗ − bdny ∗ + bmy ∗ cn − d c − x ∗ by
Subtitusi y ∗ = a b maka diperoleh C D α t V ( t ) ≤ − x ∗ − dn + m cn a − d c − x ∗ by . Jika dn + m cn < x ∗ < d c maka C D α t V ( t ) ≤ 0. Berdasarkan Lema 6 maka titik E 1 = m a + bn + d b , a b stabil asimtotik global.
Jika kondisi seperti pada Teorema 11 dipenuhi, maka untuk setiap kepadatan populasi awal bagaimanapun, kedua populasi tidak akan menuju kepunahan dan akan konstan pada saat t → ∞ . Kondisi ini yang diharapkan dalam menjaga eksistensi kedua populasi.

## 3.3. Bifurkasi Hopf
Bifurkasi Hopf merupakan suatu peristiwa perubahan kestabilan suatu titik kesetimbangan disertai munculnya limit-cycle yang mengisolasi titik kesetimbangan tersebut yang kestabilannya berlawanan dengan titik kesetimbangannya. Pada sistem dengan orde-fraksional, salah satu yang memicu terjadinya bifurkasi Hopf yaitu ketika orde dari turunannya digerakkan. Oleh karena itu, syarat terjadinya bifurkasi Hopf yaitu titik kesetimbangan tersebut memiliki sepasang nilai eigen kompleks konjugat dengan bagian real positif. Pada
model (4) , bifurkasi Hopf terjadi pada titik kesetimbangan di interior seperti yang ditunjukkan oleh teorema berikut.
Teorema 12. Jika
abm
2 ( a + bn ) 2 ! 2 < a ( bdn + ad + bm ) a + bn maka titik kesetimbangan E 1 = ( x ∗ , y ∗ ) = m a + bn + d b , a b mengalami bifurkasi Hopf pada saat α melewati α ∗ .
bukti . Tinjau kembali Teorema 10 . Karena abm
2 ( a + bn ) 2 ! 2 < a ( bdn + ad + bm ) a + bn maka matriks Jacobi dari model (4) di titik E 1 adalah sepasang nilai eigen kompleks konjugat dengan bagian real positif. Dengan mudah dapat ditunjukkan bahwa m ( α ∗ ) = 0 dan dm ( α ) d α α = α ∗ = π 2 dengan m ( α ) = απ 2 − min 1 ≤ i ≤ 2 | arg ( λ i ) | . Berdasarkan Teorema 3 di ref. [27] , E 1 mengalami bifurkasi Hopf ketika orde dari turunan α melewati titik kritis α ∗ . Titik kritis ini disebut dengan titik bifurkasi.
Dari fenomena bifurkasi Hopf, diperoleh gambaran dinamika yang lebih luas dibandingkan dengan mempelajari kestabilan titik kesetimbangan baik lokal ataupun global. Meskipun terdapat kondisi yang menunjukkan bahwa untuk nilai-nilai parameter tertentu mengakibatkan titik kesetimbangan E 1 tidak stabil, kedua populasi tidak akan menuju kepunahan, namun hanya terjadi perubahan kepadatan populasi secara periodik.

## 4. Simulasi Numerik
Simulasi numerik dilakukan dengan tujuan untuk memberikan ilustrasi dari dinamika model (4) sesuai dengan hasil analisis yang diperoleh sebelumnya. Skema yang digunakan dalam simulasi ini yaitu metode prediktor- korektor yang dikembangan oleh Diethelm dkk [28] . Dengan menggunakan metode tersebut, diperoleh skema numerik untuk model (4) sebagai berikut.
x h ( t n + 1 ) = x ( 0 ) + h α Γ ( α + 2 ) F 1 x P h ( t n + 1 ) , y P h ( t n + 1 ) + h α Γ ( α + 2 ) n ∑ j = 0 a j , n + 1 F 1 x h t j , y h t j , y h ( t n + 1 ) = y ( 0 ) + h α Γ ( α + 2 ) F 2 x P h ( t n + 1 ) , y P h ( t n + 1 ) + h α Γ ( α + 2 ) n ∑ j = 0 a j , n + 1 F 2 x h t j , y h t j ,
x P h ( t n + 1 ) = x ( 0 ) + 1 Γ ( α ) n ∑ j = 0 b j , n + 1 F 1 x h t j , y h t j y P h ( t n + 1 ) = y ( 0 ) + 1 Γ ( α ) n ∑ j = 0 b j , n + 1 F 2 x h t j , y h t j
dimana x h dan y h adalah skema korektor untuk prey dan predator , x P h dan y P h adalah skema prediktor untuk prey dan predator , h adalah step-size , dan
a j , n + 1 =    n α + 1 − ( n − α ) ( n + 1 ) α , jika j = 0, ( n − j + 2 ) α + 1 + ( n − j ) α + 1 − 2 ( n − j + 1 ) α + 1 , jika 1 ≤ j ≤ n , 1, jika j = n + 1,
b j , n + 1 = h α ( n + 1 − j ) α − ( n − j ) α .
Simulasi numerik untuk model (4) menggunakan nilai-nilai parameter yang disesuaikan dengan hasil teoritis sebelumnya yaitu sebagai berikut.
a = 0.5, b 1 = 0.3, c = 0.25, d = 0.1, m = 0.1, dan n = 0.2. (13)
Pertama-tama, dipilih nilai α = 0.92. Berdasarkan Teorema 10 dan Teorema 11 , titik kesetimbangan E 0 adalah titik pelana dan E 1 adalah stabil asimtotik. Seluruh solusi akan menuju titik kesetimbangan E 1 , lihat Gambar 1 . Seperti yang ditunjukkan pada Gambar 1b dan 1c , solusi berosilasi disekitar titik kesetimbangan E 1 dengan
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 x(t) 0 1 2 3 4 5 y( t) E 1 E 0 (a) 0 500 1000 1500 2000 2500 3000 t 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 x (t ) (b) 0 500 1000 1500 2000 2500 3000 t 0 1 2 3 4 5 y( t) (c)
Gambar 1. Simulasi numerik model (4) dengan nilai parameter pada persamaan (13) dan α = 0.92
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 x(t) 0 1 2 3 4 5 y( t) E 0 E 1 (a) 0 250 500 750 1000 1250 1500 1750 2000 t 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 x (t ) (b) 0 250 500 750 1000 1250 1500 1750 2000 t 0 1 2 3 4 5 y( t) (c)
Gambar 2. Simulasi numerik model (4) dengan nilai parameter pada persamaan (13) dan α = 0.95
amplitudo semakin mengecil dan pada akhirnya konvergen ke titik kesetimbangan E 1 . Pada saat nilai α dinaikkan menjadi 0.95, solusi disekitar E 1 yang pada awalnya stabil, menjadi tidak stabil. Berdasarkan Teorema 12 , fenomena ini perubahan kestabilan E 1 akibat perubahan nilai α disebut dengan bifurkasi Hopf yang ditunjukkan dengan munculnya sinyal periodik stabil yang mengisolasi titik kesetimbangan E 1 . Seluruh solusi yang menjauh dari titik kesetimbangan E 1 akan konvergen ke sinyal periodik tersebut. Dengan demikian, meskipun kedua titik kesetimbangan tidak stabil, tidak ada kondisi yang mengakibatkan kedua populasi punah.
Secara biologis, asumsi orde-fraksional lebih masuk akal karena sangat sulit menemukan kondisi yang ada di alam yang menyatakan bahwa jumlah dari predator dan prey selalu konstan pada periode waktu tertentu. Dengan fenomena bifurkasi tersebut, maka kepadatan dari kedua populasi akan cenderung berubah sepanjang waktu, namun tetap mempertahankan eksistensi dari kedua populasi.

## 5. Kesimpulan
Dinamika model predator-prey Lotka-Volterra orde-fraksional dengan efek Allee aditif pada predator telah dipelajari. Model ini menjamin bahwa tidak akan terjadi kepunahan terhadap kedua populasi. Meskipun populasi berkurang kepadatannya, pada akhirnya eksistensi dari kedua populasi tetap dipertahankan. Hasil analisis juga telah menunjukkan bahwa model (4) mengalami bifurkasi Hopf yang digerakkan oleh orde fraksionalnya. Solusi yang cukup dekat dengan titik kesetimbangan di interior akan menjauhi titik kesetimbangan tersebut, namun untuk waktu yang cukup panjang, kedua populasi akan berosilasi menuju ke suatu sinyal periodik. Kondisi ini menunjukkan bahwa meskipun efek Allee pada predator akan mengakibatkan perubahan kepadatan kedua populasi, eksistensi dari kedua populasi akan terjaga, baik secara konstan ataupun periodik.

## Referensi
[1] D. Indrajaya, A. Suryanto, and A. R. Alghofari, “Dynamics of modified Leslie-Gower predator-prey model with Beddington-DeAngelis functional response and additive Allee effect,” International Journal of Ecology and Development , vol. 31, no. 3, pp. 60–71, 2016.
[2] E. Rahmi and H. S. Panigoro, “Pengaruh Pemanenan terhadap Model Verhulst dengan Efek Allee,” SEMIRATA MIPAnet 2017 , no. 1, pp. 105–112, 2017.
[3] C. Arancibia-Ibarra, “The basins of attraction in a modified May–Holling–Tanner predator–prey model with Allee affect,” Nonlinear Analysis, Theory, Methods and Applications , vol. 185, pp. 15–28, 2019.
[4] L. Zhang, C. Zhang, and Z. He, “Codimension-one and codimension-two bifurcations of a discrete predator–prey system with strong Allee effect,” Mathematics and Computers in Simulation , vol. 162, pp. 155–178, 2019.
[5] W. C. Allee, Animal Aggregations, a Study in General Sociology . Chicago: University of Chicago Press, 1931.
[6] A. Suryanto, I. Darti, and S. Anam, “Stability Analysis of a Fractional Order Modified Leslie-Gower Model with Additive Allee Effect,” International Journal of Mathematics and Mathematical Sciences , vol. 2017, no. 11, pp. 1–9, 2017.
[7] T. Li, X. Huang, and X. Xie, “Stability of a stage-structured predator-prey model with allee effect and harvesting,” Communications in Mathematical Biology and Neuroscience , vol. 2019, pp. 1–20, 2019.
[8] N. Martínez-Jeraldo and P. Aguirre, “Allee effect acting on the prey species in a Leslie–Gower predation model,” Nonlinear Analysis: Real World Applications , vol. 45, pp. 895–917, 2019.
[9] H. S. Panigoro, E. Rahmi, N. Achmad, and S. L. Mahmud, “The Influence of Additive Allee Effect and Periodic Harvesting to the Dynamics of Leslie-Gower Predator-Prey Model,” Jambura Journal of Mathematics , vol. 2, no. 2, pp. 87–96, 2020.
[10] P. A. Stephens and W. J. Sutherland, “Consequences of the Allee effect for behaviour, ecology and conservation,” Trends in Ecology & Evolution , vol. 14, no. 10, pp. 401–405, oct 1999.
[11] P. Aguirre, E. Gonzalez-Olivares, and E. Saez, “Three Limit Cycles in a Leslie-Gower Predator-Prey Model with Additive Allee Effect,” SIAM J. Appl. Math. , vol. 69, no. 5, pp. 1244–1262, 2009.
[12] P. J. Pal and T. Saha, “Qualitative analysis of a predator-prey system with double Allee effect in prey,” Chaos, Solitons and Fractals , vol. 73, pp. 36–63, 2015.
[13] S. Isik, “A study of stability and bifurcation analysis in discrete-time predator-prey system involving the Allee effect,” International Journal of Biomathematics , vol. 12, no. 1, pp. 1–15, 2019.
[14] C. Rebelo and C. Soresina, “Coexistence in seasonally varying predator–prey systems with Allee effect,” Nonlinear Analysis: Real World Applications , vol. 55, p. 103140, 2020.
[15] A. J. Lotka, “Elements of Physical Biology,” Nature , vol. 116, no. 2917, pp. 461–461, 1925.
[16] A. Suryanto, I. Darti, and S. Anam, “Stability analysis of pest-predator interaction model with infectious disease in prey,” 2018, p. 020018.
[17] H. L. Li, A. Muhammadhaji, L. Zhang, and Z. Teng, “Stability analysis of a fractional-order predator–prey model
incorporating a constant prey refuge and feedback control,” Advances in Difference Equations , vol. 2018, no. 1, pp. 1–12, 2018.
[18] Z. Wang, Y. Xie, J. Lu, and Y. Li, “Stability and bifurcation of a delayed generalized fractional-order prey–predator model with interspecific competition,” Applied Mathematics and Computation , vol. 347, no. 2014, pp. 360–369, 2019.
[19] A. Suryanto, I. Darti, H. S. Panigoro, and A. Kilicman, “A Fractional-Order Predator–Prey Model with Ratio-Dependent Functional Response and Linear Harvesting,” Mathematics , vol. 7, no. 11, p. 1100, 2019.
[20] I. Petras, Fractional-order nonlinear systems: modeling, analysis and simulation . Beijing: Springer London, 2011.
[21] Z. M. Odibat and N. T. Shawagfeh, “Generalized Taylor’s formula,” Applied Mathematics and Computation , vol. 186, no. 1, pp. 286–293, 2007.
[22] D. Matignon, “Stability results for fractional differential equations with applications to control processing,” Computational engineering in systems applications , pp. 963–968, 1996.
[23] Y. Li, Y. Chen, and I. Podlubny, “Stability of fractional-order nonlinear dynamic systems: Lyapunov direct method and generalized Mittag–Leffler stability,” Computers & Mathematics with Applications , vol. 59, no. 5, pp. 1810–1821, 2010.
[24] C. Vargas-De-León, “Volterra-type Lyapunov functions for fractional-order epidemic systems,” Communications in Nonlinear Science and Numerical Simulation , vol. 24, no. 1-3, pp. 75–85, 2015.
[25] J. Huo, H. Zhao, and L. Zhu, “The effect of vaccines on backward bifurcation in a fractional order HIV model,” Nonlinear Analysis: Real World Applications , vol. 26, pp. 289–305, 2015.
[26] H. S. Panigoro, A. Suryanto, W. M. Kusumawinahyu, and I. Darti, “Dynamics of a Fractional-Order Predator-Prey Model with Infectious Diseases in Prey,” Commun. Biomath. Sci. , vol. 2, no. 2, pp. 105–117, 2019.
[27] X. Li and R. Wu, “Hopf bifurcation analysis of a new commensurate fractional-order hyperchaotic system,” Nonlinear Dynamics , vol. 78, no. 1, pp. 279–288, 2014.
[28] K. Diethelm, N. J. Ford, and A. D. Freed, “A Predictor-Corrector Approach for the Numerical Solution of Fractional Differential Equations,” Nonlinear Dynamics , vol. 29, no. 1-4, pp. 3–22, 2002.
©  2020  by  the Authors .  This  article  is  an  open  access  article  distributed  under  the  terms  and  conditions  of  the  Creative Commons  A � ribution-NonComercial  4.0  International  License. Editorial  of  JJBM: Department  of  Mathematics,  State University  of  Gorontalo,  Jln.  Prof.  Dr.  Ing.  B.  J.  Habibie,  Moutong,  Tilongkabila,  Kabupaten  Bone  Bolango,  Provinsi Gorontalo 96119, Indonesia.