MODEL PREDATOR-PREY DENGAN KONDISI INFEKSI DI KEDUA POPULASI ABSTRACT

MODEL PREDATOR-PREY DENGAN KONDISI INFEKSI DI KEDUA POPULASI Rahmazona1∗ , Khozin Mu’tamar2 1,2 Program

Studi S1 Matematika Jurusan Matematika Fakultas Matematika dan Ilmu Pengetahuan Alam Universitas Riau Kampus Bina Widya, Pekanbaru 28293 ∗ [email protected]

ABSTRACT This article discusses the analysis of mathematical models of predator-prey behavior with infectious conditions in both populations. This model uses the systems of differential equations of the model predator-prey with two compartments that are susceptible and infected in each population. The model consists of seven observed equilibrium points. Stability analysis is performed using Routh-Hurwitz criteria. Furthermore, simulation is given with specific parameters to describe the behavior and stability of the equilibrium points. Keywords: System of differential equation, model predator-prey, Routh-Hurwitz stability, infections population ABSTRAK Artikel ini membahas analisis model matematika dari perilaku predator-prey dengan kondisi infeksi di kedua populasi. Pada model ini digunakan sistem persamaan diferensial dari model predator-prey dengan dua kompartemen yaitu sehat dan terinfeksi pada masing-masing populasi. Model yang diamati terdiri atas tujuh titik ekuilibrium. Analisis kestabilan dilakukan dengan menggunakan kriteria RouthHurwitz. Selanjutnya, diberikan simulasi dengan parameter tertentu untuk menggambarkan perilaku dan kestabilan titik ekuilibriumnya. Kata kunci: Sistem persamaan diferensial, model predator-prey, kestabilan RouthHurwitz, populasi terinfeksi 1. PENDAHULUAN Dalam beberapa dekade terakhir, telah banyak penelitian tentang efek penyakit dalam sistem predator-prey. Pengaruh dari epidemi pada predasi itu pertama kali dipelajari oleh Anderson dan May [1, 2]. Mereka melakukan modifikasi persamaan Lotka-Volterra dengan pemangsa yang lebih tinggi dan tidak ada reproduksi pada 1

mangsa yang terinfeksi. Mereka menetapkan bahwa serangan penyakit cenderung untuk mengacaukan populasi mangsa-pemangsa. Hadeler dan Freedman [5] menganggap model predator-prey yang memangsa lebih mungkin pada mangsa yang terinfeksi. Mereka mempertimbangkan dalam model mereka yang pemangsa terinfeksi hanya dari mangsa yang terinfeksi. Jika dibandingkan dengan model predator-prey Lotka-Volterra, model mangsapemangsa dengan infeksi di kedua populasi lebih relevan dalam kehidupan nyata. Namun sedikit perhatian pada masalah ini, salah satu alasan utamanya adalah kurangnya mesin matematika untuk menangani meningkatnya jumlah persamaan diferensial yang terlibat dalam model predator-prey tersebut. Hal inilah yang melatarbelakangi penulis untuk mambahas lebih lanjut tentang “Model PredatorPrey dengan Kondisi Infeksi di Kedua Populasi” yang sebagian didasarkan dari artikel Bera et.al. [3]. Artikel ini tersusun atas sistematika berikut ini. Pembentukan model diberikan setelah bagian ini, dilanjutkan dengan analisis kesatabilan dan diakhiri dengan simulasi model dengan mengambil nilai parameter tertentu. 2. PEMBENTUKAN MODEL PREDATOR-PREY DENGAN KONDISI INFEKSI DI KEDUA POPULASI Dalam memodelkan predator-prey dengan kondisi infeksi di kedua populasi, populasi mangsa (prey) dibagi menjadi 2 kompartemen, yaitu populasi mangsa sehat, mangsa terinfeksi, pemangsa sehat dan pemangsa terinfeksi yang masing-masing dinotasikan S(t), I(t), X(t) dan Y (t). Beberapa asumsi yang digunakan dalam model ini yaitu (i) Penyakit hanya menyebar pada populasi yang sama yaitu penyakit mangsa hanya menyebar pada populasi mangsa saja, begitu juga dengan penyakit pemangsa hanya menyebar pada populasi pemangsa saja. (ii) Populasi yang terinfeksi tidak akan pernah pulih sehingga populasi akan mati atau tetap terinfeksi. (iii) Pemangsa yang terinfeksi tidak sanggup menangkap mangsa yang sehat atau mangsa yang sehat hanya dapat ditangkap oleh pemangsa yang sehat pula. Tetapi mangsa yang terinfeksi dapat ditangkap oleh pemangsa yang sehat maupun yang terinfeksi. (iv) Populasi pemangsa memiliki kematian alami, sedangkan populasi pemangsa yang terinfeksi penyakit mengakibatkan laju kematian yang lebih tinggi. Berdasarkan definisi dan asumsi diperoleh sistem yang menggambarkan model

2

predator-prey dengan kondisi infeksi di kedua populasi, yaitu  ( ) dS   = S 1 − S+I − ASI − BSX,  K  dt      dI  = ASI − DI − F IX − M IY,    dt dX dt dY dt

 = CSX + GIX − EXY − βX,           = EXY + N IY − (β + α)Y,  

(1)

dengan S(0) = S0 > 0, I(0) = I0 ≥ 0, X(0) = X0 > 0, Y (0) = Y0 ≥ 0. Sistem persamaan (1) memiliki dua belas parameter yang membuat analisis menjadi lebih sulit. Untuk mengurangi jumlah parameter dan untuk menentukan kombinasi variabel yang mengontrol perilaku sistem, dimisalkan s=

I X Y S , i = , x = , y = , a = AK, b = BK, c = CK, K K K K d = DK, e = EK, f = F K, g = GK, m = M K, n = N K, dan δ = β + α,

sehingga sistem persamaan (1) dapat dituliskan menjadi ds = s (1 − s − i) − asi − bsx, dt di = asi − di − f ix − miy, dt dx = csx + gix − exy − βx, dt dy = exy + niy − δy. dt

                          

(2)

3. ANALISIS KESTABILAN Dalam menganalisis kestabilan titik ekuilibrium perlu dilakukan linearisasi [8, h. 101]. Matriks linearisasi dari sistem persamaan (2) yaitu   1 − 2s − i − ai − bx −s − as −bs 0   ai as − d − f x − my −f i −mi . A=   cx gx cs + gi − ey − β −ex 0 ny ey ex + ni − δ (3) 3

Kestabilan sistem dapat dianalisis berdasarkan nilai eigen dengan menggunakan kriteria Routh-Hurwitz. Definisi 1 [6, h. 35] Misalkan nilai eigen dari matriks A adalah λi untuk i = 1, 2, . . . , n. Jika Re(λi ) < 0, maka sistem stabil asimtotik. Jika Re(λi ) ≤ 0, maka sistem stabil. Jika terdapat Re(λi ) > 0, maka sistem tidak stabil dengan Re merupakan bagian Real dari nilai eigen. Sistem persamaan (2) memiliki titik ekuilibrium [4, h. 509] berikut: (i) Titik E1 = (0, 0, 0, 0). Titik ekuilibrium pertama E1 = (0, 0, 0, 0) menghasilkan matriks A1 yaitu   1 0 0 0  0 −d 0 0   (4) A1 =   0 0 −β 0  . 0 0 0 −δ Persamaan karakteristik dari persamaan (4) diperoleh (1 − λ)(−d − λ)(−β − λ)(−δ − λ) = 0, yang dipenuhi oleh λ1 = 1 > 0, λ2 = −d < 0, λ3 = −β < 0, λ4 = −δ < 0. Berdasarkan Definisi 1, E1 tidak stabil karena terdapat λ1 = 1 > 0. (ii) Titik E2 = (1, 0, 0, 0). Selanjutnya dari persamaan (3) dengan titik ekuilibrium E2 = (1, 0, 0, 0) diperoleh bentuk matriks sebagai berikut:   −1 −(1 + a) −b 0  0 a−d 0 0  . A2 =  (5)  0 0 c−β 0  0 0 0 −δ Persamaan karakteristik dari persamaan (5) diperoleh (−1 − λ)((a − d) − λ)((c − β) − λ)(−δ − λ) = 0, yang dipenuhi oleh λ1 = −1 < 0,

λ2 = a − d,

λ3 = c − β,

λ4 = −δ < 0.

Berdasarkan Definisi 1, E2 akan stabil asimtotik jika memenuhi a < d dan c < β. 4

(iii) Titik E3 = (s∗ , 0, x∗ , 0) dengan s∗ dan x∗ diberikan oleh  β  ∗  s = , c c−β  x∗ = .  bc Kemudian, untuk titik ekuilibrium ketiga E3 = (s∗ , 0, x∗ , 0) diperoleh   bβ − βc − (1+a)β − 0 c c   β(f +ab)−c(f +bd) 0 0 0   bc A3 =  c−β g(c−β) e(c−β)  .  − b  0 − bc bc e(c−β)−bcδ 0 0 0 bc

(6)

Persamaan karakteristik dari persamaan (6) adalah (λbc − ec + eβ + bcδ)(λ3 + B1 λ2 + B2 λ + B3 ) = 0, dengan

c(f + bd) − β(f + ab) + bβ , bc −β 2 (f + ab) + cβ(f + bd) + bc2 β − bcβ 2 B2 = , bc2 β 3 (f + ab) + c2 β(f + bd) − cβ 2 (2f + ab + bd) B3 = , bc2

B1 =

sehingga diperoleh nilai eigen λ1 =

e(c − β) − bcδ < 0, bc

jika eβ + bcδ > ec.

Berdasarkan Kriteria Routh-Hurwitz [7, h. 222], E3 akan stabil asimtotik jika memenuhi B1 > 0, B3 > 0 dan B1 B2 > B3 . (iv) Titik E4 = (s∗ , 0, x∗ , y ∗ ) dengan s∗ , x∗ , y ∗ diberikan oleh  e − bδ  ∗  s = ,   e  δ ∗ x = ,  e   ce − bcδ − eβ  ∗  . y = 2 e Titik ekuilibrium keempat E4 = (s∗ , 0, x∗ , y ∗ ) menghasilkan matriks linearisasi A4 yaitu   −s∗ −(1 + a)s∗ −bs∗ 0  0 as∗ − d + f x∗ − my ∗ 0 0  . (7) A4 =  ∗  cx∗ gx 0 −ex∗  0 ny ∗ ey ∗ 0 5

Persamaan karakteristik dari persamaan (7) diperoleh λ4 + C1 λ3 + C2 λ2 + C3 λ + C4 = 0, dengan C1 C2 C3 C4

= s∗ − as∗ − f x∗ + my ∗ + d, = e2 x∗ y ∗ − s∗ (s∗ − C1 ), = bcs∗ x∗ (C1 − s∗ ) + e2 x∗ y ∗ C1 , = e2 s∗ x∗ y ∗ (C1 − s∗ ).

Berdasarkan Kriteria Routh-Hurwitz [7, h. 222], E4 akan stabil asimtotik jika memenuhi C1 > 0, C3 > 0, C4 > 0 dan C1 C2 C3 > C32 + C12 C4 . (v) Titik E5 = (s∗ , i∗ , 0, 0) dengan s∗ dan i∗ diberikan oleh  d  ∗  s = , a a−d  i∗ = .  a(1 + a) Selanjutnya titik ekuilibrium kelima E5 = (s∗ , i∗ , 0, 0) diperoleh bentuk matriks sebagai berikut:   −s∗ −(1 + a)s∗ −bs∗ 0  ai∗ 0 −f i∗ −mi∗  . A5 =  (8) ∗ ∗  0  0 cs + gi − β 0 0 0 0 ni∗ − δ Persamaan karakteristik dari persamaan (8) diperoleh (−λ + ni∗ − δ)(−λ + cs∗ − gi∗ − β)(λ2 + s∗ λ + ai∗ s∗ (1 + a)) = 0, sehingga n(a − d) − δ, a(1 + a) (1 + a)(cd − aβ) − g(a − d) λ2 = , a(1 + a) √ −d ± (1 − 4a)d2 − 4a2 d λ3,4 = . 2a λ1 =

Berdasarkan Definisi 1, titik ekuilibrium E5 akan stabil asimtotik jika memenuhi n(a − d) < aδ(1 + a) dan (1 + a)(cd − aβ) < g(a − d).

6

(vi) Titik E6 = (s∗ , i∗ , x∗ , 0) dengan s∗ , i∗ , x∗ diberikan oleh s∗ i∗ x∗

g(f + bd) − β(f + af ) , = g(f + ab) − c(f + af ) β(f + ab) − c(f + bd) = , g(f + ab) − c(f + af ) ag(f + bd) − dg(f + ab) + (f + af )(cd − aβ) = . gf (f + ab) − cf (f + af )

              

Kemudian, untuk titik ekuilibrium keenam E6 = (s∗ , i∗ , x∗ , 0) diperoleh   −s∗ −(1 + a)s∗ −bs∗ 0   ai∗ as∗ − d − f x∗ −f i∗ −mi∗ . A6 =  ∗ ∗ ∗ ∗   cx∗ gx cs + gi − β −ex 0 0 0 ex∗ + ni∗ − δ

(9)

Persamaan karakteristik dari persamaan (9) adalah (λ − ex∗ − ni∗ + δ)(λ3 + D1 λ2 + D2 λ + D3 ) = 0, dengan D1 = s∗ (1 − a − c) − gi∗ + f x∗ + d + β, D2 = s∗ x∗ (c (b − f ) + f ) + i∗ s∗ (g (a − 1) + a (a + 1)) + (s∗ )2 (c (a − 1) − a) + s∗ (d(1 − c) + β(1 + a)) + β(f x∗ + d) − dgi∗ , D3 = s∗ (x∗ (bcd + βf ) − s∗ (cd + βa) + i∗ x∗ (abg − cf (a + 1)) + (i∗ )2 ag(a − 1) + i∗ (βa(a + 1) − dg) + ai∗ s∗ (g − c(a + 1)) − s∗ x∗ (cf + abc) + bcf ∗ (x∗ )2 + ac(s∗ )2 + βd), sehingga diperoleh nilai eigen λ1 = ex∗ + ni∗ − δ < 0,

jika ex∗ + ni∗ < δ.

Berdasarkan Kriteria Routh-Hurwitz [7, h. 222], E6 akan stabil asimtotik jika memenuhi D1 > 0, D3 > 0 dan D1 D2 > D3 . (vii) Titik E7 = (s∗ , i∗ , x∗ , y ∗ ) dengan s∗ , i∗ , x∗ , y ∗ diberikan oleh P2 − (1 + a)P1 − bP3 s = , P2 P1 i∗ = , P2 P3 x∗ = , P2 (g − c − ac)P1 + (c − β)P2 − bcP3 y∗ = , eP2 ∗

                    

7

dan P1 = e(abδ + f δ + cm + de) − (bcmδ + emβ + ae2 ), P2 = e(abn + f n + (1 + δ)cm+) − (bcmn + emg + (1 + a)ae2 , P3 = (aen + mnβ + (1 + a)cmδ) − (cmn + den + (1 + a)aeδ + gmδ). Terakhir, untuk titik ekuilibrium ketujuh E7  −s∗ −(1 + a)s∗  ai∗ as∗ − d − f x∗ − my ∗ A7 =   cx∗ gx∗ 0 ny ∗

= (s∗ , i∗ , x∗ , y ∗ ) diperoleh  −bs∗ 0  −f i∗ −mi∗ . ∗  0 −ex ∗ ∗ ∗ ey ex + ni − δ

(10)

Persamaan karakteristik dari persamaan (10) adalah λ4 + F1 λ3 + F2 λ2 + F3 λ + F4 = 0, dengan F1 = s∗ (1 − a) − ni∗ − x∗ (e − f ) + my ∗ + d + δ, F2 = δ(f x∗ + my ∗ + s∗ (1 − a) + d) − ni∗ (s∗ (1 − a) + d + f x∗ ) + ex∗ (s∗ (a − 1) + f x∗ + y ∗ (e − m) − d) + s∗ (a(i∗ + ai∗ − s∗ ) + x∗ (f + bc) + my ∗ + d) + gf i∗ x∗ , F3 = bs∗ c(x∗ )2 (f − e) + f i∗ s∗ x∗ (n + g − c(1 + a)) − es∗ x∗ (a(ey ∗ − s∗ + i∗ (1 + ai∗ )) + y ∗ (m − e) + f x∗ + d ) + bcs∗ x∗ ( δ − as∗ + d − ni∗ + my ∗ ) + ex∗ y ∗ ( gmi∗ + de + f ex∗ + mey ∗ ) + s∗ ( δ ( f x∗ + my ∗ − as∗ + d ) − dni∗ ) + gx∗ f i∗ (i∗ n + δ), F4 = s∗ x∗ (f i∗ ex∗ (c(1 + a) − g) + f i∗ (ac(ni∗ − δ) − n(ey ∗ + gi∗ + c(bx∗ − i∗ )) + δ(g − c)) + ei∗ y ∗ (m(g − c) + a (e + ai∗ − bn − cm)) − ex∗ (b ( agi∗ + cmy ∗ + cd + cf x∗ ) − ef y ∗ ) + δb ( agi∗ + c( my ∗ + f x∗ − as∗ + d )) + bni∗ (agi∗ + c − acs∗ (as∗ − d)) + ey ∗ (e(my ∗ + d) − as∗ y ∗ )). Berdasarkan Kriteria Routh-Hurwitz [7, h. 222], E7 akan stabil asimtotik jika memenuhi F1 > 0, F3 > 0, F4 > 0 dan F1 F2 F3 > F32 + F12 F4 .

4. SIMULASI MODEL Misalkan dipilih parameter dengan nilai a= 0.4, b = 1, c = 0.1, d = 0.5, e = 0.2, f = 0.2, g = 0.15, m = 1.5, n = 1.4, δ = 0.5, β = 0.2 dan misalkan diberikan nilai awal s(0) = 0.5, i(0) = 0.5, x(0) = 0.5 dan y(0) = 0.5. Dengan menggunakan MAPLE hanya diperoleh dua titik ekuilibrium yang bernilai positif yaitu E1 = (0, 0, 0, 0) dan E2 = (1, 0, 0, 0) sedangkan untuk titik ekuilibrium lainnya tidak memenuhi syarat. Selanjutnya dilakukan simulasi model dengan menggunakan dua titik ekuilibrium tersebut yaitu E1 = (0, 0, 0, 0) dan E2 = (1, 0, 0, 0). 8

Titik Ekuilibrium pertama E1 = (0, 0, 0, 0) Perhatikan kembali matriks linearisasi untuk titik ekuilibrium pertama yaitu E1 = (0, 0, 0, 0) pada persamaan (4), sehingga diperoleh   1 0 0 0  0 −0.5 0 0  , A1 =   0 0 −0.2 0  0 0 0 −0.5 sehingga diperoleh nilai eigen sebagai berikut: λ1 λ2 λ3 λ4

 = 1,    = −0.2, = −0.5,    = −0.5.

(11)

Solusi khusus untuk titik ekuilibrium pertama sebagai berikut:  s(t) = 0.5et ,    i(t) = 0.5e−0.5t , x(t) = 0.5e−0.2t ,    y(t) = 0.5e−0.5t .

(12)

Sistem persamaan (2) memiliki nilai eigen pada persamaan (11). Sistem persamaan (2) di sekitar E1 = (0, 0, 0, 0) bersifat tidak stabil, sehingga sistem akan menjauhi titik E1 dari nilai awal yang diberikan seperti yang ditunjukkan pada Gambar 1. i 0.5 0.4 0.3 0.2 0.1 0

0.2

0.4

0.6

0.8 s

Gambar 1: Perilaku sistem di sekitar E1 . Solusi khusus sistem persamaan (2) di titik ekuilibrium E1 seperti pada persamaan (12) ditunjukkan pada Gambar 2. Pada Gambar 2 terlihat bahwa untuk t yang sangat besar, solusi khusus sistem persamaan (2) di titik E1 akan menuju (∞, 0, 0, 0) sehingga sistem memiliki solusi yang tidak stabil. Karena untuk t menuju tak hingga solusi sistem persamaan (2) menuju 9

200 150 100 50

0

2

4

6t

Gambar 2: Grafik solusi sistem di titik E1 . (∞, 0, 0, 0), maka jumlah populasi mangsa sehat akan terus bertambah karena nilai s yang tidak stabil menuju tak hingga sehingga jumlah populasi mangsa terinfeksi terus berkurang karena adanya kapasitas daya tampung yang ada. Perbandingan jumlah populasi tersebut ditunjukkan pada Gambar 3. i 0.5 0.4

0.3

0.2 0

2

4

6 s

Gambar 3: Perbandingan grafik s dan i pada E1 . Titik Ekuilibrium Kedua E2 = (1, 0, 0, 0) Perhatikan kembali matriks linearisasi untuk titik ekuilibrium E2 pada persamaan (5). Karena titik ekuilibrium E2 tak nol, maka dilakukan pergeseran seperti pada persamaan (13) berikut: x′ = A(x − xe ). Misalkan θ = x − xe , θe = 0 dan x′ = θ ′ diperoleh θ ′ linearisasi dari E2 diperoleh sebagai berikut:   ′   θs −1 −1.4 −1 0   θi′   0 −0.1 0 0    ′ =  θx   0 0 −0.1 0  ′ θy 0 0 0 −0.5

(13) = Aθ sehingga matriks  θs θi  , θx  θy

10

dengan 

 −1 −1.4 −1 0  0 −0.1 0 0  , A2 =   0 0 −0.1 0  0 0 0 −0.5 sehingga diperoleh nilai eigen berikut:  λ1 = −1,    λ2 = −0.5, λ3 = −0.1,    λ4 = −0.1.

(14)

Solusi khusus E2 diperoleh sebagai berikut: θs = 1.83333e−t − 1.33333e−0.1t , θi = 0.5e−0.1t , θx = 0.5e−0.1t , θy = 0.5e−0.5t .

      

Karena θ = x − xe , maka x = θ + xe sehingga diperoleh solusi khusus untuk titik ekuilibrium E2 = (1, 0, 0, 0) sebagai berikut:  s(t) = 1 + 1.83333e−t − 1.33333e−0.1t ,    i(t) = 0.5e−0.1t , (15) x(t) = 0.5e−0.1t ,    y(t) = 0.5e−0.5t . Sistem persamaan (2) memiliki nilai eigen pada persamaan (14). Sistem persamaan (2) di sekitar E2 = (1, 0, 0, 0) bersifat stabil asimtotik, sehingga sistem akan mendekati titik E2 dari nilai awal yang diberikan seperti yang ditunjukkan pada Gambar 4. i 0.5 0.4 0.3 0.2 0.1 0

0.5

1 s

Gambar 4: Perilaku sistem di sekitar E2 . Solusi khusus sistem persamaan (2) di titik ekuilibrium E2 seperti pada persamaan (15) ditunjukkan pada Gambar 5. Pada Gambar 5 terlihat bahwa untuk 11

t yang sangat besar maka solusi khusus sistem persamaan (2) di titik E2 akan menuju (1, 0, 0, 0) sistem memiliki solusi yang stabil asimtotik. 1

0.5

0

10

20

30

40

50 t

Gambar 5: Grafik solusi sistem di titik E2 . Karena untuk t menuju tak hingga solusi sistem persamaan (2) menuju (1, 0, 0, 0), maka jumlah populasi mangsa sehat akan terus bertambah karena nilai s yang menuju suatu nilai sebanding dengan jumlah populasi mangsa terinfeksi juga terus bertambah. Perbandingan jumlah populasi tersebut ditunjukkan pada Gambar 6. i 0.5 0.49 0.48 0.47 0.46 0.2

0.3

0.4

0.5 s

Gambar 6: Perbandingan grafik s dan i pada E2 .

DAFTAR PUSTAKA [1] R. M. Anderson dan R. M. May, Population biology of infectious diseases, PartI, Nature, 280 (1979), 361-367. [2] R. M. Anderson dan R. M. May, The invasion, persistence, and spread of infectious diseases within animal and plant communities, Philosophical Transactions of the Royal Society of London, B314 (1986), 533-570. [3] S. P. Bera, A. Maiti dan G. P. Samanta, A prey-predator model with infection in both prey and predator, Filomat, 29 (2015), 1753-1767.

12

[4] W. E. Boyce dan R. C. DiPrima, Elementary Differential Equation and Boundary Value Problem, Tenth Edition, Jonh Wiley and Sons, New York, 2012. [5] K. P. Hadeler dan H. I. Freedman, Predator-prey populations with parasitic infection, Journal of Mathematical Biology, 27 (1989), 609-631. [6] X. Liao, L. Wang dan P. Yu, Stability of Dynamical System, Elsevier, Amsterdam, 2007. [7] K. Ogata, Modern Control Engineering, Fifth Edition, Pearson, New Jersey, 2010. [8] L. Perko, Differential Equations and Dynamical Systems, Second Edition, Springer-Verlag, New York, 1996.

13