Limits: Journal of Mathematics and Its Applications E-ISSN: 2579-8936 P-ISSN: 1829-605X Vol. 17, No. 2, Desember 2020, 157-172 DOI: http://dx.doi.org/10.12962/limits.v17i2.6901

## Analisis Dinamik pada Model Kanker Serviks dengan Vaksinasi dan Screening
Karunia Theda Kristanti 1 *, Trisilowati 2 , Agus Widodo 3
1,2,3 Fakultas Matematika dan Ilmu Pengetahuan Alam, Universitas Brawijaya, Malang, Indonesia e-mail: kristantikarunia@gmail.com
Diajukan: 11 Mei 2020,  Diperbaiki: 19 Juli 2020, Diterima:27 Juli 2020

## Abstrak
Pada  paper  ini  dibahas  analisis  dinamik  model  penyebaran  kanker  serviks  dengan  melibatkan tindakan  vaksinasi  dan screening .  Penyebab  utama  terjadinya  kanker  serviks  adalah  karena  seseorang terinfeksi Human Papillomavirus (HPV). Infeksi ini dapat menular karena adanya kontak langsung melalui hubungan seksual antara subpopulasi wanita rentan dengan pria terinfeksi HPV maupun kontak langsung antara  pria  rentan  dengan  wanita  terinfeksi  HPV.  Pada  model  ini  diasumsikan  vaksin  diberikan  pada subpopulasi  wanita  rentan  saja  dengan  salah  satu  jenis  vaksin.  Sementara  itu, screening dilakukan  oleh subpopulasi wanita terifeksi HPV sebagai upaya deteksi dini untuk mencegah terjadinya kanker serviks. Hasil  analisis  dinamik  menunjukkan  bahwa  model  penyebaran  kanker  serviks  dengan  vaksinasi  dan screening memiliki  dua  titik  kesetimbangan  yaitu  titik  kesetimbangan  bebas  penyakit  dan  titik kesetimbangan endemi. Eksistensi dan kestabilan lokal titik kesetimbangan bergantung pada nilai angka reproduksi  dasar 𝑅 0 .  Berdasarkan  hasil  analisis,  titik  kesetimbangan  bebas  penyakit  eksis  tanpa  syarat, sedangkan titik kesetimbangan endemi eksis jika 𝑅 0 > 1. Titik kesetimbangan bebas penyakit bersifat stabil asimtotik  lokal  jika  nilai 𝑅 0 < 1 dan  titik  kesetimbangan  endemi  bersifat  stabil  asimtotik  lokal  jika memenuhi  kriteria  Routh-Hurwitz.  Simulasi  numerik  yang  dilakukan  mendukung  hasil  analisis  dinamik yang diperoleh.
Kata Kunci : analisis kestabilan lokal, angka reproduksi dasar, Human Papillomavirus , kanker serviks.

## Abstract
In this paper, the dynamical analysis of cervical cancer model with vaccination and screening is discussed. The leading cause of cervical cancer is because of Human Papillomavirus (HPV). This infection can  be  transmitted  due  to  direct  contact  through  sexual  intercourse  between  susceptible  females  with infected males and  contact directly between  susceptible males with infected  females. In this model,  it is assumed  that  vaccination  is  given  to  susceptible  females  with  any  of  various  vaccines.  Meanwhile, screening is done by infected females with HPV as an early detection to prevent cervical cancer. Dynamical analysis results show that the model has two equilibrium points, namely the disease-free equilibrium point and  endemic  equilibrium  point.  The  existence  and  stability  of  equilibrium  points  depend  on  the  basic reproduction  number 𝑅 0 . The  disease-free  equilibrium  point  exists  unconditionally,  while  the  endemic equilibrium  exists  when 𝑅 0 > 1 .  If 𝑅 0 < 1 ,  then  the  disease-free  equilibrium  point  is  local asymptomatically stable. The endemic equilibrium point is locally asymptomatically stable if it satisfies the Routh-Hurwitz criteria. The numerical simulation is illustrated to support the analysis result.
Keywords : local stability analysis, basic reproduction number, Human Papillomavirus, cervical cancer .

## 1 Pendahuluan
Kanker serviks adalah kanker leher rahim yang ditandai dengan perubahan sel abnormal dan secara bertahap berkembang menjadi kanker. Penyakit ini dapat dicegah jika terdeteksi lebih awal dan dilakukan tindakan yang tepat, tetapi menjadi penyakit yang mematikan jika terlambat dalam penangannya  [1].  Kanker  serviks  terkonfirmasi  70%  disebabkan  oleh  infeksi  Human Papillomavirus  (HPV).  Selain  itu,  penyebab  kanker  serviks  didukung  karena  hubungan  seksual yang  dilakukan  terlalu  dini,  usia  persalinan  pertama  yang  belum  cukup  umur,  rendahnya pendidikan, dan jumlah pasangan seksual yang berganti-ganti [2].
Vaksinasi dapat dipilih sebagai salah satu cara untuk merangsang sistem kekebalan tubuh dalam memproduksi antibodi sebagai upaya pencegahan terhadap HPV yang dapat menginfeksi sel  dalam  tubuh.  Terdapat  tiga jenis  vaksin  yang  dapat  digunakan  untuk  mencegah  penyebaran HPV yaitu vaksin bivalen, vaksin quadrivalen, dan vaksin nonavalen [3], ketiga jenis vaksisn ini dapat digunakan secara bersamaan [4]. Pemberian vaksin secara lengkap dibagi dalam tiga kali dosis dalam kurun waktu 6 bulan. Jika pemberian vaksin tidak diberikan dengan lengkap dapat memengaruhi efikasi vaksin dalam mencegah infeksi HPV [5]. Sementara itu, screening bertujuan untuk mencari adanya masalah kesehatan yang mana jika tidak dilakukan dengan benar menjadi tidak efektif. Screening dapat dilakukan dengan tes Papsmear atau IVA. Tes Papsmear merupakan pemeriksaan  sitologi  pada  serviks  guna  mendeteksi  adanya  perubahan  pada  daerah  serviks, sedangkan  IVA  merupakan  tes  visual  menggunakan  larutan  asam  cuka  dan  iosium  lugol  pada serviks dan melihat perubahan warna yang terjadi untuk melihat adanya sel yang displasi [6].
Model matematika telah banyak digunakan untuk mendeskripsikan dampak dari vaksinasi dan screening terhadap  infeksi  HPV  yang  dapat  berakibat  kanker  serviks.  Pada  tahun  2006, Elbasha  telah  menganalisis  dampak  vaksin quadrivalen terhadap  penyebaran  HPV  dengan mengembangkan  model  yang  membagi  populasi  menjadi  dua  yaitu  pria  dan  wanita  yang didasarkan  pada  model SIR  (Susceptible-Infected-Recovery) [7] . Penelitian  tentang  penyebaran HPV dengan vaksinasi semakin berkembang. Majed dan Lounes telah menganalisis penyebaran HPV  dengan  model SIS  (Susceptible-Infected-Susceptible) dengan  menerapkan  vaksin quadrivalen [8]. Sementara itu, Van de Velde dkk., membandingkan keefektifan vaksin bivalen , quadrivalen , dan nonavalen terhadap upaya pencegahan penyebaran infeksi HPV [9]. Selanjutnya, Malik dkk. (2016) menerapkan kontrol optimal untuk memperoleh strategi pengendalian infeksi HPV  dengan  ketiga  jenis  vaksin,  dimana  vaksinasi  hanya  diberikan  pada  populasi  wanita  [4]. Selain vaksinasi, penyebaran HPV dengan melibatkan screening telah dibahas oleh Saldana dkk., dengan tujuan untuk mengurangi jumlah populasi terinfeksi HPV [10].

## Karunia Theda Kristanti, Trisilowati, Agus Widodo
Model matematika penyebaran HPV dengan melibatkan kanker serviks telah dibahas oleh beberapa  peneliti.  Hasnawati  pada  tahun  2017  membentuk  model  penyebaran  kanker  serviks dengan  melibatkan  subpopulasi  vaksinasi  dan  subpopulasi  wanita  terkena  kanker  serviks  serta subpopulasi wanita sembuh bisa kembali menjadi wanita rentan [11]. Berbeda dengan model milik Puspita 2019 menunjukkan adanya interaksi antara subpopulasi wanita rentan dengan subpopulasi wanita terinfeksi tanpa adanya vaksinasi [12]. Pada model milik puspitasari, subpopulasi wanita terinfeksi dibedakan menjadi subpopulasi wanita terkena kanker serviks dan tidak terkena kanker serviks.  Sementara  itu,  ditahun  yang  sama  Sado  membentuk  model  penyebaran  kanker  serviks dengan  menambahkan  kompartemen  subpopulasi  wanita  tervaksin.  Model  penyebaran  kanker serviks  pada  penelitian  ini  melibatkan  kompartemen  wanita  dan  pria  dengan  melibatkan subpopulasi  penderita  kanker  serviks.  Selain  itu,  vaksinasi dan screening dipilih  sebagai  upaya untuk mencegah penyebaran infeksi HPV untuk menekan terjadinya kanker serviks.

## 2 Metode Penelitian
Pada  tahap  awal  penelitian  dilakukan  konstruksi  model  kanker  serviks  dengan  adanya tindakan  vaksinasi  dan screening .  Selanjutnya  menentukan  titik  kesetimbangan  model  yang diperoleh pada saat sistem berada pada keadaan setimbang dengan menentukan syarat eksistensi dari  titik  kesetimbangan  tersebut.  Pada  tahap  berikutnya  dilakukan  analisis  kestabilan  titik kesetimbangan  model.  Penentuan  kestabilan  diawali  dengan  linearisasi  model  yang  selanjutnya ditentukan akar-akar persamaan karakterikstik pada sistem linear hampiran yang telah dibentuk. Kestabilan  titk  kesetimbangan  ditentukan  menggunakan  kriteria  Routh-Hourwitz.  Pada  tahap akhir  dilakukan  simulasi  numerik  untuk  mengilustrasikan  hasil  analisis  dinamik  dengan menggunakan perangkat lunak MATLAB.

## 3 Hasil dan Pembahasan
Pada bagian ini dibahas konstruksi model, penentuan titik kesetimbangan, angka reproduksi dasar, analisis kestabilan lokal, serta simulasi numerik untuk mengilustrasikan hasil analisis yang telah diperoleh.

## 3.1 Konstruksi Model
Model  matematika  yang  dibentuk  pada  penelitian  ini  merupakan  penyebaran  HPV  yang dapat  menyebabkan  kanker  serviks  pada  wanita.  Model  terdiri  dari  kompartemen  wanita  dan kompartemen  pria.  Infeksi  HPV  terjadi  karena  adanya  hubungan  seksual  atau  interaksi  antara subpopulasi  wanita  rentan  dengan  subpopulasi  pria  terinfeksi  HPV  demikian  pula  sebaliknya.
Vaksinasi diberikan  pada  subpopulasi  wanita rentan, tetapi hal  ini  tidak  menutup  kemungkinan wanita yang sudah divaksin dapat terinfeksi HPV. Wanita aktif secara seksual memiliki peluang yang  tinggi  terinfeksi  HPV,  oleh  karena  itu screening dapat  dilakukan  untuk  mengetahui  ada tidaknya  sel  abnormal  yang  memiliki  peluang  kanker  serviks.  Sementara  itu,  wanita  terinfeksi HPV yang tidak melakukan screening akan memiliki peluang lebih tinggi terkena kanker serviks karena tidak adanya perlindungan maupun pengobatan yang dilakukan. Infeksi HPV yang terjadi baik pada subpopulasi wanita dan pria dapat disembuhkan dengan tindakan pengobatan. Berikut adalah diagram kompartemen yang menggambarkan penyebaran HPV pada pria dan wanita yang aktif secara seksual sampai terjadinya kanker serviks pada subpopulasi wanita.
Gambar 1. Diagram kompartemen model penyebaran kanker serviks
Berdasarkan  diagram  kompartemen  pada  Gambar  1,  model  penyebaran  kanker  serviks  dapat dituliskan dalam bentuk persamaan diferensial sebagai berikut.
𝑑𝑆 𝑓 (𝑡)
𝑑𝑡 = 𝜋 𝑓 𝑁 𝑓 (𝑡) − 𝜙𝑆 𝑓 (𝑡) − 𝑐 𝑓 𝛽 𝑚 𝐽(𝑡)𝑆 𝑓 (𝑡) 𝑁 𝑚 (𝑡) − 𝜇 𝑓 𝑆 𝑓 (𝑡),
𝑑𝑉(𝑡) 𝑑𝑡 =  𝜙𝑆 𝑓 (𝑡) − 𝜇 𝑓 𝑉(𝑡) − (1 − 𝜀) 𝑐 𝑓 𝛽 𝑚 𝐽(𝑡)𝑉(𝑡) 𝑁 𝑚 , 𝑑𝐼 ℎ (𝑡) 𝑑𝑡 = 𝑐 𝑓 𝛽 𝑚 𝐽(𝑡)𝑆 𝑓 (𝑡) 𝑁 𝑚 (𝑡) + (1 − 𝜀) 𝑐 𝑓 𝛽 𝑚 𝐽(𝑡)𝑉(𝑡) 𝑁 𝑚 (𝑡) − 𝑝𝐼 ℎ (𝑡) − 𝜇 𝑓 𝐼 ℎ (𝑡) − 𝜃(1 − 𝑝)𝐼 ℎ (𝑡),
𝑑𝐼 𝑠𝑐 (𝑡) 𝑑𝑡 = 𝑝𝐼 ℎ (𝑡) − 𝜎𝐼 𝑠𝑐 (𝑡) − 𝜇 𝑓 𝐼 𝑠𝑐 (𝑡) − 𝛿𝐼 𝑠𝑐 (𝑡),
𝑑𝐶(𝑡) 𝑑𝑡 = 𝛿𝐼 𝑠𝑐 (𝑡) + θ(1 − 𝑝)𝐼 ℎ (𝑡) − (𝜇 𝑓 + 𝛾)𝐶(𝑡),
𝑑𝑅 𝑓 (𝑡)
𝑑𝑡 = 𝜎𝐼 𝑠𝑐 (𝑡) − 𝜇 𝑓 𝑅 𝑓 (𝑡),
𝑑𝑆 𝑚 (𝑡) 𝑑𝑡 = 𝜋 𝑚 𝑁 𝑚 (𝑡) − 𝑐 𝑚 𝛽 𝑓 𝐼 ℎ (𝑡)𝑆 𝑚 (𝑡) 𝑁 𝑓 (𝑡) − 𝜇 𝑚 𝑆 𝑚 (𝑡) + 𝛼𝑅 𝑚 (𝑡), (1)
𝑑𝐽(𝑡) 𝑑𝑡 = 𝑐 𝑚 𝛽 𝑓 𝐼 ℎ (𝑡)𝑆 𝑚 (𝑡) 𝑁 𝑓 (𝑡) − 𝜁𝐽(𝑡) − 𝜇 𝑚 𝐽(𝑡),
𝑑𝑅 𝑚 (𝑡) 𝑑𝑡 = 𝜁𝐽(𝑡) − 𝜇 𝑚 𝑅 𝑚 (𝑡),
dengan
𝑁 𝑓 (𝑡) = 𝑆 𝑓 (𝑡) + 𝑉(𝑡) + 𝐼 ℎ (𝑡) + 𝐼 𝑠𝑐 (𝑡) + 𝐶(𝑡) + 𝑅 𝑓 (𝑡), dan 𝑁 𝑚 (𝑡) = 𝑆 𝑚 (𝑡) + 𝐽(𝑡) + 𝑅 𝑚 (𝑡).
Parameter yang digunakan pada model diuraikan pada Tabel 1 sebagai berikut.
Tabel 1. Nilai Parameter Parameter Deskripsi 𝜋 𝑓 (𝜋 𝑚 ) Laju populasi wanita (pria) aktif secara seksual
𝜇 𝑓 (𝜇 𝑚 )
Kematian alami wanita (pria)
𝛽 𝑓 (𝛽 𝑚 ) Peluang infeksi berpindah dari wanita ke pria (pria ke wanita)
𝑐 𝑓 (𝑐 𝑚 ) Rata-rata banyaknya kontak seksual 𝜙 Tingkat vaksinasi 𝜀 Efikasi Vaksin 𝑝 Laju Screening 𝜃 Laju wanita terinfeksi HPV terkena kanker serviks 𝜎 Tingkat populasi wanita sembuh 𝜁 Tingkat populasi pria sembuh 𝛿 Laju subpopulasi wanita screening yang terkena kanker
serviks
𝛾 Laju kematian karena kanker serviks
Sistem persamaan (1) dapat disederhanakan dengan melakukan penskalaan. Misalkan,
𝑆̅ 𝑓 (𝑡) = 𝑆 𝑓 (𝑡) 𝑁 𝑓 (𝑡) , 𝑉̅(𝑡) = 𝑉(𝑡) 𝑁 𝑓 (𝑡) , 𝐼̅ ℎ (𝑡) = 𝐼 ℎ (𝑡) 𝑁 𝑓 (𝑡) , 𝐼̅ 𝑠𝑐 (𝑡) = 𝐼 𝑠𝑐 (𝑡) 𝑁 𝑓 (𝑡) , 𝐶̅(𝑡) = 𝐶(𝑡) 𝑁 𝑓 (𝑡) , 𝑅̅ 𝑓 (𝑡) = 𝑅 𝑓 (𝑡) 𝑁 𝑓 (𝑡) , 𝑆̅ 𝑚 (𝑡) = 𝑆 𝑚 (𝑡) 𝑁 𝑚 (𝑡) , 𝐽̅(𝑡) = 𝐽(𝑡) 𝑁 𝑚 (𝑡) , 𝑅̅ 𝑚 (𝑡) = 𝑅 𝑚 (𝑡) 𝑁 𝑚 (𝑡) .
Hasil penskalaan di atas sebagai berikut.
𝑑𝑆̅ 𝑓 (𝑡) 𝑑𝑡 = 𝜋 𝑓 − 𝜙𝑆̅ 𝑓 (𝑡) − 𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑆̅ 𝑓 (𝑡) − 𝜋 𝑓 𝑆̅ 𝑓 (𝑡),
𝑑𝑉̅(𝑡) 𝑑𝑡 = 𝜙𝑆̅ 𝑓 (𝑡) − (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑉̅(𝑡) − 𝜋 𝑓 𝑉̅(𝑡),
𝑑𝐼̅ ℎ (𝑡) 𝑑𝑡 = 𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑆̅ 𝑓 (𝑡) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑉̅(𝑡) − 𝑝𝐼̅ ℎ (𝑡) − 𝜃(1 − 𝑝)𝐼̅ ℎ (𝑡) − 𝜋 𝑓 𝐼̅ ℎ (𝑡),
𝑑𝐼̅ 𝑠𝑐 (𝑡) 𝑑𝑡 = 𝑝𝐼̅ ℎ − 𝜎𝐼̅ 𝑠𝑐 (𝑡) − 𝛿𝐼̅ 𝑠𝑐 (𝑡) − 𝜋 𝑓 𝐼̅ 𝑠𝑐 (𝑡),
𝑑𝐶̅(𝑡) 𝑑𝑡 = 𝛿𝐼̅ 𝑠𝑐 (𝑡) + θ(1 − 𝑝)𝐼̅ ℎ (𝑡) − 𝛾𝐶̅(𝑡) − 𝜋 𝑓 𝐶̅(𝑡),
𝑑𝑅̅ 𝑓 (𝑡)
𝑑𝑡 = 𝜎𝐼̅ 𝑠𝑐 (𝑡) − 𝜋 𝑓 𝑅̅ 𝑓 (𝑡),
(2)
𝑑𝑆̅ 𝑚 (𝑡) 𝑑𝑡 = 𝜋 𝑚 − 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡)𝑆̅ 𝑚 (𝑡) − 𝜋 𝑚 𝑆̅ 𝑚 (𝑡),
𝑑𝐽̅(𝑡) 𝑑𝑡 = 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡)𝑆̅ 𝑚 (𝑡) − 𝜁𝐽̅(𝑡) − 𝜋 𝑚 𝐽̅(𝑡),
𝑑𝑅̅ 𝑚 (𝑡) 𝑑𝑡 = 𝜁𝐽̅(𝑡) − 𝜋 𝑚 𝑅̅ 𝑚 (𝑡).
Berdasarkan sistem (2) tampak bahwa 𝑅̅ 𝑓 (𝑡) dan 𝑅̅ 𝑚 (𝑡) tidak terlibat di dalam persamaan yang lain, sehingga sistem (2) dapat direduksi untuk memudahkan dalam analisis. Berikut adalah sistem persamaan hasil reduksi.
𝑑𝑆̅ 𝑓 (𝑡)
𝑑𝑡 = 𝜋 𝑓 − 𝜙𝑆̅ 𝑓 (𝑡) − 𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑆̅ 𝑓 (𝑡) − 𝜋 𝑓 𝑆̅ 𝑓 (𝑡),
𝑑𝑉̅(𝑡)
𝑑𝑡 = 𝜙𝑆̅ 𝑓 (𝑡) − (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑉̅(𝑡) − 𝜋 𝑓 𝑉̅(𝑡),
𝑑𝐼̅ ℎ (𝑡) 𝑑𝑡 = 𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑆̅ 𝑓 (𝑡) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑉̅(𝑡) − 𝑝𝐼̅ ℎ (𝑡) − 𝜃(1 − 𝑝)𝐼̅ ℎ (𝑡) − 𝜋 𝑓 𝐼̅ ℎ (𝑡),
𝑑𝐼̅ 𝑠𝑐 (𝑡) 𝑑𝑡 = 𝑝𝐼̅ ℎ − 𝜎𝐼̅ 𝑠𝑐 (𝑡) − 𝛿𝐼̅ 𝑠𝑐 (𝑡) − 𝜋 𝑓 𝐼̅ 𝑠𝑐 (𝑡),
𝑑𝐶̅(𝑡) 𝑑𝑡 = 𝛿𝐼̅ 𝑠𝑐 (𝑡) + θ(1 − 𝑝)𝐼̅ ℎ (𝑡) − 𝛾𝐶̅(𝑡) − 𝜋 𝑓 𝐶̅(𝑡),
𝑑𝑆̅ 𝑚 (𝑡) 𝑑𝑡 = 𝜋 𝑚 − 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡)𝑆̅ 𝑚 (𝑡) − 𝜋 𝑚 𝑆̅ 𝑚 (𝑡),
𝑑𝐽̅(𝑡) 𝑑𝑡 = 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡)𝑆̅ 𝑚 (𝑡) − 𝜁𝐽̅(𝑡) − 𝜋 𝑚 𝐽̅(𝑡).
3.2  Titik Kesetimbangan dan Angka Reproduksi Dasar
Titik kesetimbangan sistem (3) diperoleh ketika
𝑑𝑆̅ 𝑓 (𝑡) 𝑑𝑡 = 𝑑𝑉̅(𝑡) 𝑑𝑡 = 𝑑𝐼̅ ℎ (𝑡) 𝑑𝑡 =, 𝑑𝐼̅ 𝑠𝑐 (𝑡) 𝑑𝑡 =, 𝑑𝐶̅(𝑡) 𝑑𝑡 =, 𝑑𝑆̅ 𝑚 (𝑡) 𝑑𝑡 = 𝑑𝐽̅(𝑡) 𝑑𝑡 = 0.
Berdasarkan perhitungan analisis diperoleh dua titik kesetimbangan yaitu,
1. Titik kesetimbangan bebas penyakit
𝐸 0 = (𝑆̅ 𝑓 (𝑡) 0 , 𝑉̅(𝑡) 0 , 𝐼̅ ℎ (𝑡) 0 , 𝐼̅ 𝑠𝑐 (𝑡) 0 , 𝐶̅(𝑡) 0 , 𝑆̅ 𝑚 (𝑡) 0 , 𝐽̅(𝑡) 0 ) = ( 𝜋 𝑓 𝜙 + 𝜋 𝑓 , 𝜙 𝜙 + 𝜋 𝑓 , 0,0,0,1,0).
Selanjutnya  ditentukan  angka  reproduksi  dasar (𝑅 0 ) yang  menjadi  angka  penting  untuk mengetahui terjadi tidaknya penyebaran penyakit pada suatu populasi. 𝑅 0 menyatakan rata-rata banyaknya individu yang terinfeksi akibat satu individu yang telah terinfeksi sebelumnya. Matriks Generasi Selanjutnya dapat digunakan untuk menentukan 𝑅 0 [13].
(3)
Misalkan 𝑥 = (𝐼̅ ℎ , 𝐼̅ 𝑠𝑐 , 𝐶̅, 𝐽̅) ,  maka  berdasarkan  model  pada  persamaan  (3)  dapat  dibentuk matriks 𝐹 dan 𝑉 yang terdiri dari subpopulasi terinfeksi sebagai berikut.
𝑑𝑥 𝑑𝑡 = 𝐹(𝑥) − 𝑉(𝑥),
dengan,
𝐹(𝑥) = ( 𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑆̅ 𝑓 (𝑡) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡)𝑉̅(𝑡)
0 0 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡)𝑆̅ 𝑚 (𝑡) ),
dan
V(x) = ( 𝑉 1 𝑉 2 𝑉 3 𝑉 4 ) = ( (𝑝 + 𝜋 𝑓 + 𝜃(1 − 𝑝)) 𝐼̅ ℎ (𝑡)
−𝑝𝐼̅ ℎ (𝑡) + (𝜎 + 𝛿 + 𝜋 𝑓 )𝐼̅ 𝑠𝑐 (𝑡) −𝛿𝐼̅ 𝑠𝑐 (𝑡) − θ(1 − 𝑝)𝐼̅ ℎ (𝑡) + (𝛾 + 𝜋 𝑓 )𝐶̅(𝑡) (𝜁 + 𝜋 𝑚 )𝐽̅(𝑡) ) .
Matriks Jacobi untuk 𝐹 dan 𝑉 pada titik kesetimbangan bebas penyakit 𝐸 0 adalah sebagai berikut.
𝐷𝐹(𝐸 0 ) =
( 0 0 0 𝑐 𝑓 𝛽 𝑚 𝜙 + 𝜋 𝑓 (𝜋 𝑓 + (1 − 𝜀)𝜙) 0 0 0 0 0 0 0 0 𝑐 𝑚 𝛽 𝑓 0 0
0 ) ,
𝐷𝑉(𝐸 0 ) = ( 𝑝 + 𝜋 𝑓 + 𝜃(1 − 𝑝) 0 0 0 −𝑝 𝜎 + 𝛿 + 𝜋 𝑓 0 0 −𝜃(1 − 𝑝) −𝛿 𝛾 + 𝜋 𝑓 0 0 0 0 𝜁 + 𝜋 𝑚 ) .
Selanjutnya  angka  reproduksi  dasar (𝑅 0 ) diperoleh  dari  radius spectral matriks (𝐷𝐹(𝐸 0 ))(𝐷𝑉 −1 (𝐸 0 )) atau nilai eigen terbesar dari matriks (𝐷𝐹(𝐸 0 ))(𝐷𝑉 −1 (𝐸 0 )) , yaitu
𝑅 0 = √ 𝑐 𝑚 𝛽 𝑓 𝑐 𝑓 𝛽 𝑚 (𝜋 𝑓 + (1 − 𝜀)𝜙) (𝑝 + 𝜋 𝑓 + 𝜃(1 − 𝑝))(𝜙 + 𝜋 𝑓 )(𝜁 + 𝜋 𝑚 ) .

## 2. Titik kesetimbangan endemi
Titik kesetimbangan endemi ( 𝐸 ∗ ) diperoleh ketika 𝐼̅ ℎ (𝑡) ≠ 0 sehingga,
𝐸 ∗ = (𝑆̅ 𝑓 (𝑡) ∗ , 𝑉̅(𝑡) ∗ , 𝐼̅ ℎ (𝑡) ∗ , 𝐼̅ 𝑠𝑐 (𝑡) ∗ , 𝐶̅(𝑡) ∗ , 𝑆̅ 𝑚 (𝑡) ∗ , 𝐽̅(𝑡) ∗ ),
dengan
𝑆̅ 𝑓 (𝑡) ∗ = 𝜋 𝑓 𝜋 𝑚 (𝜁 + 𝜋 𝑚 ) + 𝜋 𝑓 𝑐 𝑚 𝛽 𝑓 (𝜁 + 𝜋 𝑚 )𝐼̅ ℎ (𝑡) ∗ (𝜙 + 𝜋 𝑓 )(𝜁𝜋 𝑚 + 𝜋 𝑚 2 ) + ((𝜙 + 𝜋 𝑓 )(𝜁𝑐 𝑚 𝛽 𝑓 + 𝜋 𝑚 𝑐 𝑚 𝛽 𝑓 ) + 𝑐 𝑓 𝛽 𝑚 𝜋 𝑚 𝑐 𝑚 𝛽 𝑓 ) 𝐼̅ ℎ (𝑡) ∗ ,
164
𝑉̅(𝑡) ∗ = 𝜙𝜋 𝑓 (𝜁𝜋 𝑚 + 𝜋 𝑚 2 + (𝜁𝑐 𝑚 𝛽 𝑓 + 𝜋 𝑚 𝑐 𝑚 𝛽 𝑓 )𝐼̅ ℎ (𝑡) ∗ ) 2 (𝜋 𝑓 𝜋 𝑚 2 + 𝜋 𝑓 𝜋 𝑚 𝜁 + 𝜓 1 𝐼̅ ℎ (𝑡)) (𝜋 𝑚 2 𝜙 + 𝜋 𝑚 𝜙𝜁 + 𝜋 𝑓 𝜋 𝑚 𝜁 + 𝜋 𝑓 𝜋 𝑚 2 + 𝜓 2 𝐼̅ ℎ (𝑡) ∗ ) , 𝐼̅ 𝑠𝑐 (𝑡) ∗ = 𝑝 (𝜎 + 𝜋 𝑓 + 𝛿) 𝐼̅ ℎ (𝑡) ∗ ,
𝐶̅(𝑡) ∗ = [ 𝛿𝑝 + θ(1 − 𝑝)(𝜎 + 𝜋 𝑓 + 𝛿) (𝛾 + 𝜋 𝑓 )(𝜎 + 𝜋 𝑓 + 𝛿) ] 𝐼̅ ℎ (𝑡) ∗ ,
𝑆̅ 𝑚 (𝑡) ∗ = 𝜋 𝑚 𝜋 𝑚 + 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡) ∗ ,
𝐽̅(𝑡) ∗ = 𝜋 𝑚 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡) ∗ (𝜁 + 𝜋 𝑚 )(𝜋 𝑚 + 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡) ∗ ) .
𝐼̅ ℎ (𝑡) ∗ = −𝑞 2 + √𝑞 2 2 − 4𝑞 1 𝑞 3 2𝑞 1 ,
dengan,
𝑞 1 = (𝑐 𝑚 𝛽 𝑓 ) 2 [𝜋 𝑓 𝜁 + 𝜋 𝑓 𝜋 𝑚 + 𝜙𝜁 + 𝜙𝜋 𝑚 + 𝑐 𝑓 𝛽 𝑚 𝜋 𝑚 ][𝜋 𝑓 𝜁 + 𝜋 𝑓 𝜋 𝑚 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝜋 𝑚 ] [𝑝 + 𝜋 𝑓 + 𝜃(1 − 𝑝)] > 0,
𝑞 2 = 𝜋 𝑚 𝑐 𝑚 𝛽 𝑓 (𝜋 𝑓 + 𝜙)(𝜁 + 𝜋 𝑚 )[𝜋 𝑓 (𝜁 + 𝜋 𝑚 ) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝜋 𝑚 ][𝑝 + 𝜋 𝑓 + 𝜃(1 − 𝑝)]
+𝑐 𝑚 𝛽 𝑓 𝜋 𝑓 𝜋 𝑚 (𝜁 + 𝜋 𝑚 )[(𝜋 𝑓 + 𝜙)(𝜁 + 𝜋 𝑚 ) + 𝑐 𝑓 𝛽 𝑚 𝜋 𝑚 ][𝑝 + 𝜋 𝑓 + 𝜃(1 − 𝑝)]
−𝜋 𝑓 𝑐 𝑓 𝛽 𝑚 𝜋 𝑚 𝑐 𝑚 2 𝛽 𝑓 2 [𝜋 𝑓 (𝜁 + 𝜋 𝑚 ) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝜋 𝑚 + (1 − 𝜀)𝜙(𝜁 + 𝜋 𝑚 )],
𝑞 3 = 𝜋 𝑚 2 𝜋 𝑓 (𝜁 + 𝜋 𝑚 ) 2 (𝜋 𝑓 + 𝜙) (𝑝 + 𝜋 𝑓 + 𝜃(1 − 𝑝)) [1 − 𝑅 0 2 ].
Eksistensi titik kesetimbangan endemi (𝐸 ∗ ) bergantung pada nilai 𝑅 0 . Jika 𝑅 0 > 1 maka diperoleh nilai 𝑞 3 < 0 sehingga  diperoleh  solusi  tunggal  untuk 𝐼̅ ℎ (𝑡) ∗ yang  bernilai  positif.  Hal  ini mengakibatkan 𝑆̅ 𝑓 (𝑡) ∗ , 𝑉̅(𝑡) ∗ , 𝐼̅ 𝑠𝑐 (𝑡) ∗ , 𝐶̅(𝑡) ∗ , 𝑆̅ 𝑚 (𝑡) ∗ , 𝐽̅(𝑡) ∗ semuanya  bernilai  positif  dan  titik kesetimbangan endemi eksis.
3.3 Kestabilan lokal titik kesetimbangan

## 3.3.1  Kestabilan Titik Kesetimbangan E 0
Kestabilan  titik  kesetimbangan  bebas  penyakit 𝐸 0 diperoleh  dengan  menentukan  nilai eigen dari Matriks Jacobi sistem (3) ketika 𝐸 0 . Titik kesetimbangan 𝐸 0 bersifat stabil asimtotik lokal jika bagian real dari nilai eigen bernilai negatif.
Matriks Jacobi di 𝐸 0 = ( 𝜋 𝑓 𝜙+𝜋 𝑓 , 𝜙 𝜙+𝜋 𝑓 , 0,0,0,1,0) adalah
𝒥(𝐸 0 ) = [ −(𝜙 + 𝜋 𝑓 ) 0 0 0 0 0 − 𝑐 𝑓 𝛽 𝑚 𝜋 𝑓 𝜙 + 𝜋 𝑓 𝜙 −𝜋 𝑓 0 0 0 0 − (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝜙 𝜙 + 𝜋 𝑓 0 0 −(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) 0 0 0 𝑐 𝑓 𝛽 𝑚 𝜙 + 𝜋 𝑓 (𝜋 𝑓 + (1 − 𝜀)𝜙) 0 0 𝑝 −(𝜎 + 𝛿 + 𝜋 𝑓 ) 0 0 0 0 0 𝜃(1 − 𝑝) 𝛿 −(𝛾 + 𝜋 𝑓 ) 0 0 0 0 −𝑐 𝑚 𝛽 𝑓 0 0 −𝜋 𝑚 0 0 0 𝑐 𝑚 𝛽 𝑓 0 0 0 −(𝜁 + 𝜋 𝑚 ) ]
Nilai  eigen  dari  matriks 𝒥(𝐸 0 ) dapat  diperoleh  dengan  cara  menyelesaikan |𝒥(𝐸 0 ) − 𝑟𝐼| = 0, sehingga diperoleh nilai eigen yaitu 𝑟 1 = −𝜋 𝑓 , 𝑟 2 = −(𝜙 + 𝜋 𝑓 ), 𝑟 3 = −𝜋 𝑚 , 𝑟 4 = −(𝛾 + 𝜋 𝑓 ), 𝑟 5 = −(𝜎 + 𝛿 + 𝜋 𝑓 ) sementara 𝑟 6 dan 𝑟 7 ditentukan dengan menyelesaikan trace (𝐿) dan |𝐿| , dimana
𝐿 = ( −(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) 𝑐 𝑓 𝛽 𝑚 𝜙 + 𝜋 𝑓 (𝜋 𝑓 + (1 − 𝜀)𝜙) 𝑐 𝑚 𝛽 𝑓 −(𝜁 + 𝜋 𝑚 ) ) = 0.
Dengan menggunakan metode trace dan determinan [14], diperoleh nilai 𝑡𝑟(𝐿) < 0 dan |𝐿| > 0 maka  didapatkan 𝑟 6 < 0 dan 𝑟 7 < 0 jika √ 𝑐 𝑚 𝛽 𝑓 𝑐 𝑓 𝛽 𝑚 (𝜋 𝑓 +(1−𝜀)𝜙) (𝑝+𝜋 𝑓 +𝜃(1−𝑝))(𝜙+𝜋 𝑓 )(𝜁+𝜋 𝑚 ) < 1 atau 𝑅 0 < 1. Oleh karena itu dapat disimpulkan bahwa titik kesetimbangan bebas penyakit 𝐸 0 bersifat stabil asimtotik lokal jika 𝑅 0 < 1.
3.3.2  Kestabilan Titik Kesetimbangan 𝐸 ∗
Kestabilan  titik  kesetimbangan  bebas  penyakit 𝐸 ∗ diperoleh  dengan  menentukan  nilai eigen  dari  matriks  jacobi  sistem  (3)  ketika 𝐸 ∗ . Titik  kesetimbangan 𝐸 ∗ bersifat  stabil  asimtotik lokal jika memenuhi kriteria Routh-Hurwitz [15]. Matriks Jacobi di 𝐸 ∗ = (𝑆̅ 𝑓 (𝑡) ∗ , 𝑉̅(𝑡) ∗ , 𝐼̅ ℎ (𝑡) ∗ , 𝐼̅ 𝑠𝑐 (𝑡) ∗ , 𝐶̅(𝑡) ∗ , 𝑆̅ 𝑚 (𝑡) ∗ , 𝐽̅(𝑡) ∗ ), adalah
𝒥(𝐸 ∗ ) = [ 𝑗 11 0 0 0 0 0 𝑗 17 𝑗 21 𝑗 22 0 0 0 0 𝑗 27 𝑗 31 𝑗 32 𝑗 33 0 0 0 𝑗 37 0 0 𝑗 43 𝑗 44 0 0 0 0 0 𝑗 53 𝑗 54 𝑗 55 0 0 0 0 𝑗 63 0 0 𝑗 66 0 0 0 𝑗 73 0 0 𝑗 76 𝑗 77 ] ,
dengan,
𝑗 11 = −(𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡) ∗ + 𝜋 𝑓 ),          𝑗 37 = 𝑐 𝑓 𝛽 𝑚 (𝑆̅ 𝑓 (𝑡) ∗ + (1 − 𝜀)𝑉̅(𝑡) ∗ ), 𝑗 73 = 𝑐 𝑚 𝛽 𝑓 𝑆̅ 𝑚 (𝑡) ∗ , 𝑗 17 = −𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 (𝑡) ∗ ,                               𝑗 43 = 𝑝,                                                             𝑗 76 = 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡) ∗ , 𝑗 21 = 𝜙,                                                    𝑗 44 = −(𝜎 + 𝛿 + 𝜋 𝑓 ),                                   𝑗 77 = −(𝜁 + 𝜋 𝑚 ). 𝑗 22 = − ((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡) ∗ + 𝜋 𝑓 ),   𝑗 53 = 𝜃(1 − 𝑝),
𝑗 27 = −(1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅(𝑡) ∗ ,                  𝑗 54 = 𝛿,
𝑗 31 = 𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡) ∗ ,                                    𝑗 55 = −(𝛾 + 𝜋 𝑓 ),
𝑗 32 = (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅(𝑡) ∗ ,                      𝑗 63 = −𝑐 𝑚 𝛽 𝑓 𝑆̅ 𝑚 (𝑡) ∗ ,
𝑗 33 = −(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ),            𝑗 66 = −(𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ (𝑡) ∗ + 𝜋 𝑚 ),
Nilai  eigen  dari  matriks 𝒥(𝐸 ∗ ) dapat  diperoleh  dengan  cara  menyelesaikan |𝒥(𝐸 ∗ ) − 𝑟𝐼| = 0, sehingga diperoleh nilai eigen yaitu 𝑟 1 = −(𝛾 + 𝜋 𝑓 ), 𝑟 2 = −(𝜎 + 𝛿 + 𝜋 𝑓 ), dan |𝐺| = 0,
dengan
|𝐺| = | | 𝑗 11 − 𝑟 0 0 0 𝑗 17 𝑗 21 𝑗 22 − 𝑟 0 0 𝑗 27 𝑗 31 𝑗 32 𝑗 33 − 𝑟 0 𝑗 37 0 0 𝑗 63 𝑗 66 − 𝑟 0 0 0 𝑗 73 𝑗 76 𝑗 77 − 𝑟 | | ,
memenuhi  kriteria  Routh-Hurwitz.  Dari  perhitungan |𝐺| = 0 maka  didapatkan  persamaan karakteristik sebagai berikut.
𝑟 5 + 𝑎 1 𝑟 4 + 𝑎 2 𝑟 3 + 𝑎 3 𝑟 2 + 𝑎 4 𝑟 + 𝑎 5 = 0 ,
dengan
𝑎 1 = 2𝜋 𝑚 + 3𝜋 𝑓 + 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜁 + 𝜙 + 𝑝 + 𝜃(1 − 𝑝) + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡),
𝑎 2 = (𝜁 + 𝜋 𝑚 + 𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 3𝜋 𝑓 + 𝑝 + 𝜃(1 − 𝑝) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡)) (𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 ) + (𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 3𝜋 𝑓 + 𝑝 + 𝜃(1 − 𝑝) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡)) (𝜁 + 𝜋 𝑚 )
+(𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) (𝑝 + 𝜃(1 − 𝑝) + 2𝜋 𝑓 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡))
+(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) ((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) + (𝑐 𝑚 𝛽 𝑓 𝑆̅ 𝑚 ∗ (𝑡)) (𝑐 𝑓 𝛽 𝑚 (𝑆̅ 𝑓 ∗ (𝑡) + (1 − 𝜀)𝑉̅ ∗ (𝑡))),
𝑎 3 = (𝜁 + 𝜋 𝑚 )(𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 ) (𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 3𝜋 𝑓 + 𝑝 + 𝜃(1 − 𝑝) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡))
+(𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 )(𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 ) (𝑝 + 𝜃(1 − 𝑝) + 2𝜋 𝑓 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡))
+(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) ((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) (𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 ) +(𝜁 + 𝜋 𝑚 )(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) (𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 2𝜋 𝑓 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡)) +(𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) ((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) (𝜁 + 𝜋 𝑚 + 𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 )
+ (𝑐 𝑚 𝛽 𝑓 𝑆̅ 𝑚 ∗ (𝑡)) (2𝜋 𝑓 𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡) + (1 − 𝜀)𝑐 𝑓 2 𝛽 𝑚 2 𝐽̅ ∗ (𝑡) (𝑆̅ 𝑓 ∗ (𝑡) + 𝑉̅ ∗ (𝑡)) +𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡)𝜙 + 2𝜋 𝑓 (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡) + (1 − 𝜀)𝜙𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡)) +𝑐 𝑚 𝛽 𝑓 𝑆̅ 𝑚 ∗ (𝑡)𝜋 𝑚 𝑐 𝑓 𝛽 𝑚 (𝑆̅ 𝑓 ∗ (𝑡) + (1 − 𝜀)𝑉̅ ∗ (𝑡)),
(4)
𝑎 4 = (𝜁𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜁𝜋 𝑚 + 𝜋 𝑚 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 2 )((𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 )
(𝑝 + 𝜃(1 − 𝑝) + 2𝜋 𝑓 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡)) + (𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 )
((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) + (𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 )(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) ((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) (𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 2𝜋 𝑚 + 𝜁) + 𝑐 𝑚 𝛽 𝑓 𝑆̅ 𝑚 ∗ (𝑡)(2𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡)𝜋 𝑚 𝜋 𝑓 +(1 − 𝜀)𝑐 𝑓 2 𝛽 𝑚 2 𝑆̅ 𝑓 ∗ (𝑡)𝐽̅ ∗ (𝑡)𝜋 𝑚 + 𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡)𝜋 𝑚 𝜙 +2(1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡)𝜋 𝑚 𝜋 𝑓 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡)𝜋 𝑚 𝜙 +(1 − 𝜀)𝑐 𝑓 2 𝛽 𝑚 2 𝑉̅ ∗ (𝑡)𝜋 𝑚 𝐽̅ ∗ (𝑡) + 𝜙𝜋 𝑓 𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡)𝜙𝜋 𝑓
+(1 − 𝜀)𝜋 𝑓 𝑐 𝑓 2 𝛽 𝑚 2 𝐽̅ ∗ (𝑡)𝑉̅ ∗ (𝑡) + (1 − 𝜀) 2 𝑐 𝑓 3 𝛽 𝑚 3 𝐽̅ ∗ (𝑡) 2 𝑉̅ ∗ (𝑡) +𝜋 𝑓 2 𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡) + 𝜋 𝑓 2 (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡) + 𝜋 𝑓 (1 − 𝜀)𝑐 𝑓 2 𝛽 𝑚 2 𝐽̅ ∗ (𝑡)𝑆̅ 𝑓 ∗ (𝑡)).
𝑎 5 = 𝜋 𝑚 𝑐 𝑚 𝛽 𝑓 𝑆̅ 𝑚 ∗ (𝑡)(𝜙𝜋 𝑓 𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡) + 𝜙𝜋 𝑓 (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡) + 𝜋 𝑓 (1 − 𝜀)𝑐 𝑓 2 𝛽 𝑚 2 𝐽̅ ∗ (𝑡)𝑉̅ ∗ (𝑡)
+𝜋 𝑓 2 𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡) + 𝜋 𝑓 2 (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡) + 𝜋 𝑓 (1 − 𝜀)𝑐 𝑓 2 𝛽 𝑚 2 𝐽̅ ∗ (𝑡)𝑆̅ 𝑓 ∗ (𝑡))
+(𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 )(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) ((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 )
(𝜁 + 𝜋 𝑚 )(𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 ),
Berdasarkan kriteria Routh Hurwitz akar-akar persamaan karakteristik (4) memiliki bagian real negatif karena memenuhi kriteria sebagai berikut.
1. 𝑎 1 = 2𝜋 𝑚 + 3𝜋 𝑓 + 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜁 + 𝜙 + 𝑝 + 𝜃(1 − 𝑝) + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) > 0, 𝑎 5 = 𝜋 𝑚 𝑐 𝑚 𝛽 𝑓 𝑆̅ 𝑚 ∗ (𝑡)(𝜙𝜋 𝑓 𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡) + 𝜙𝜋 𝑓 (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡) + 𝜋 𝑓 (1 − 𝜀)𝑐 𝑓 2 𝛽 𝑚 2 𝐽̅ ∗ (𝑡)𝑉̅ ∗ (𝑡)
+𝜋 𝑓 2 𝑐 𝑓 𝛽 𝑚 𝑆̅ 𝑓 ∗ (𝑡) + 𝜋 𝑓 2 (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝑉̅ ∗ (𝑡) + 𝜋 𝑓 (1 − 𝜀)𝑐 𝑓 2 𝛽 𝑚 2 𝐽̅ ∗ (𝑡)𝑆̅ 𝑓 ∗ (𝑡))
+(𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 )(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) ((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) (𝜁 + 𝜋 𝑚 )(𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 ) > 0,
2. 𝑎 1 𝑎 2 − 𝑎 3 = (𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 ) 2 (𝜁 + 𝜋 𝑚 + 𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 3𝜋 𝑓 + 𝑝 + 𝜃(1 − 𝑝) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡))
+2(𝜁 + 𝜋 𝑚 )(𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 ) (𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 3𝜋 𝑓 + 𝑝 + 𝜃(1 − 𝑝) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡))
+2(𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 )(𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 ) (𝑝 + 𝜃(1 − 𝑝) + 2𝜋 𝑓 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡))
+2(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) ((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) (𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 2𝜋 𝑚 + 𝜁 + 𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) +(𝜁 + 𝜋 𝑚 ) 2 (𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 𝜋 𝑚 + 𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 3𝜋 𝑓 + 𝑝 + 𝜃(1 − 𝑝) + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡))
+𝑐 𝑚 𝛽 𝑓 𝑆̅ 𝑚 ∗ (𝑡) (𝑐 𝑓 𝛽 𝑚 (𝑆̅ 𝑓 ∗ (𝑡) + (1 − 𝜀)𝑉̅ ∗ (𝑡))) (𝜁 + 𝜋 𝑚 + 𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 + 𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡)) +(𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) 2 (𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 2𝜋 𝑚 + 𝜁 + 𝑝 + 𝜃(1 − 𝑝) + 2𝜋 𝑓 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡))
+2(𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 )(𝜁 + 𝜋 𝑚 ) (𝑝 + 𝜃(1 − 𝑝) + 2𝜋 𝑓 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡))
+(𝑝 + 𝜃(1 − 𝑝) + 𝜋 𝑓 ) 2 (𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 2𝜋 𝑚 + 𝜁 + 𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 2𝜋 𝑓 + (1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡)) + ((1 − 𝜀)𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 𝜋 𝑓 ) 2 (𝑐 𝑚 𝛽 𝑓 𝐼̅ ℎ ∗ (𝑡) + 2𝜋 𝑚 + 𝜁 + 𝜙 + 𝑐 𝑓 𝛽 𝑚 𝐽̅ ∗ (𝑡) + 2𝜋 𝑓 + 𝑝 + 𝜃(1 − 𝑝)) +𝑐 𝑚 𝛽 𝑓 𝑐 𝑓 2 𝛽 𝑚 2 𝑆̅ 𝑚 ∗ (𝑡)𝐽̅ ∗ (𝑡) (𝑆̅ 𝑓 ∗ (𝑡) + (1 − 𝜀) 2 𝑉̅ ∗ (𝑡)) > 0 3. 𝑎 1 𝑎 2 𝑎 3 + 𝑎 1 𝑎 5 − 𝑎 3 2 − 𝑎 1 2 𝑎 4 > 0,
4. 𝑎 1 𝑎 2 𝑎 3 𝑎 4 + 2𝑎 1 𝑎 4 𝑎 5 + 𝑎 2 𝑎 3 𝑎 5 − 𝑎 1 2 𝑎 4 2 − 𝑎 1 𝑎 2 2 𝑎 5 − 𝑎 3 2 𝑎 4 − 𝑎 5 2 > 0.

## Analisis Dinamik pada Model Kanker Serviks dengan Vaksinasi dan Screening
Kriteria  (1)  dan  (2)  telah  terpenuhi  dengan  detail  perhitungan  yang  telah  dijabarkan.  Dengan demikian, titik kesetimbangan endemi 𝐸 ∗ stabil asimtotik lokal jika memenuhi kriteria (3) dan (4) yang dibuktikan secara numerik.

## 3.4 Simulasi Numerik
Pada  bagian  ini  dilakukan  simulasi  numerik  untuk  mengilustrasikan  kestabilan  titik kesetimbangan  bebas  penyakit 𝐸 0 ketika 𝑅 0 < 1 dan  juga  dilakukan  simulasi  numerik  untuk mengilustrasikan  kestabilan  titik  kesetimbangan  endemi 𝐸 ∗ ketika 𝑅 0 > 1 .  Simulasi  numerik dilakukan untuk mendukung hasil analisis yang telah diperoleh. Nilai awal yang digunakan adalah 𝑆̅ 𝑓 (0) = 0.3, 𝑉̅(0) = 0.2, 𝐼̅ ℎ (0) = 0.1, 𝐼̅ 𝑠𝑐 (0) = 0.1, 𝐶̅(0) = 0.2,  𝑆̅ 𝑚 (0) = 0.5, 𝐽̅(0) = 0.3, dan nilai parameter yang digunakan pada simulasi ini tercantum pada Tabel 2.
Tabel 2. Nilai parameter untuk simulasi numerik Parameter Nilai Parameter Simulasi 𝑅 0 < 1 Nilai Parameter Simulasi 𝑅 0 > 1 𝜋 𝑓 (𝜋 𝑚 ) 0.5 0.5
𝛽 𝑓 (𝛽 𝑚 ) 0.4 0.8 𝑐 𝑓 (𝑐 𝑚 ) 2 3 𝜙 0.05 0.05 𝜀 0.95 0.95 𝑝 0.01 0.01 𝜃 0.1 0.1 𝜎 0.7 0.7 𝜁 0.7 0.7 𝛿 0.05 0.05 𝛾 0.01 0.01 3.4.1 Simulasi I: Kestabilan titik 𝐸 0 saat 𝑅 0 < 1
Pada simulasi ini digunakan nilai parameter seperti pada Tabel 2 dengan 𝛽 𝑓 = 0.4, 𝛽 𝑚 = 0.4, 𝑐 𝑓 = 2, dan 𝑐 𝑚 = 2, sehingga didapatkan nilai 𝑅 0 = 0.8945 < 1 .
Gambar 2. Solusi numerik kestabilan 𝐸 0 saat 𝑅 0 = 0.8945 < 1

## Karunia Theda Kristanti, Trisilowati, Agus Widodo
Berdasarkan Gambar 2 ditunjukkan bahwa solusi numerik konvergen menuju titik kesetimbangan 𝐸 0 (𝑆̅ 𝑓 , 𝑉̅, 𝐼̅ ℎ , 𝐼̅ 𝑠𝑐 , 𝐶̅, 𝑆̅ 𝑚 , 𝐽̅) = (0.9091,0.0909,0,0,0,1,0) .  Hasil  simulasi  tersebut  sesuai  dengan hasil  analisis  yang  diperoleh  bahwa  jika 𝑅 0 < 1 tidak  ada  infeksi  pada  subpopulasi  dan  titik kesetimbangan bebas penyakit 𝐸 0 bersifat stabil asimtotik.
3.4.2 Simulasi I: Kestabilan titik 𝐸 ∗ saat 𝑅 0 > 1
Pada simulasi ini digunakan nilai parameter seperti pada Tabel 2 dengan 𝛽 𝑓 = 0.8, 𝛽 𝑚 = 0.8, 𝑐 𝑓 = 3, dan 𝑐 𝑚 = 3, sehingga didapatkan nilai 𝑅 0 = 2.6835 > 1 .
Gambar 3. Solusi numerik kestabilan 𝐸 ∗ saat 𝑅 0 = 2.6835 < 1
Berdasarkan Gambar 3 ditunjukkan bahwa solusi numerik konvergen menuju titik kesetimbangan 𝐸 ∗ (𝑆̅ 𝑓 , 𝑉̅, 𝐼̅ ℎ , 𝐼̅ 𝑠𝑐 , 𝐶̅, 𝑆̅ 𝑚 , 𝐽̅) = (0.4041, 0.0516, 0.4583, 0.0037, 0.0893, 0.3125,0.2864) . Sementara itu  juga  diperoleh  nilai 𝑎 1 𝑎 2 𝑎 3 + 𝑎 1 𝑎 5 − 𝑎 3 2 − 𝑎 1 2 𝑎 4 = 337.7289 > 0 dan 𝑎 1 𝑎 2 𝑎 3 𝑎 4 + 2𝑎 1 𝑎 4 𝑎 5 + 𝑎 2 𝑎 3 𝑎 5 − 𝑎 1 2 𝑎 4 2 − 𝑎 1 𝑎 2 2 𝑎 5 − 𝑎 3 2 𝑎 4 − 𝑎 5 2 = 1481.8 > 0 sehingga kriteria Routh- Hurwitz terpenuhi. Dengan  demikian titik kesetimbangan 𝐸 ∗ bersifat stabil asimtotik lokal.  Oleh karena  itu  dapat  disimpulkan  bahwa  hasil  simulasi  sesuai  dengan  hasil  analisis  yang  diperoleh bahwa  jika 𝑅 0 > 1 maka  terjadi  penyebaran  infeksi  pada  subpopulasi  dan  titik  kesetimbangan endemi 𝐸 ∗ bersifat stabil asimtotik.
Gambar 4. Grafik solusi jumlah subpopulasi terinfeksi saat 𝑝 = 0, 𝜙 = 0, 𝑝 = 0.01, 𝜙 = 0.05
Gambar  4  menunjukkan  pengaruh  tidak  dan  adanya  pemberian  vaksinasi  maupun screening terhadap jumlah subpopulasi wanita terinfeksi ( 𝐼̅ ℎ ), subpopulasi wanita terkena kanker serviks ( 𝐶̅ ), dan subpopulasi pria terinfeksi ( 𝐽̅ ). Parameter 𝑝 = 0, 𝜙 = 0 menunjukkan tidak ada tindakan  vaksinasi  dan screening yang  diberikan.  Pada  Gambar  4  terlihat  bahwa  pemberian vaksinasi  dan screening dengan  parameter 𝑝 = 0.01, 𝜙 = 0.05 mampu  mengurangi  jumlah subpopulasi  terinfeksi  diantaranya  mengurangi  jumlah  subpopulasi  wanita  terinfeksi  HPV, subpopulasi wanita terkena kanker serviks, dan subpopulasi pria terinfeksi HPV.

## 4 Simpulan
Berdasarkan hasil pembahasan pada penelitian ini, maka dapat disimpulkan sebagai berikut.
 Model  penyebaran  kanker  serviks  dengan  vaksinasi  dan screening mempunyai  dua  titik kesetimbangan. Titik kesetimbangan tersebut yaitu titik kesetimbangan bebas penyakit dan titik kesetimbangan endemi. Eksistensi dan kestabilan titik kesetimbangan bergantung pada nilai angka reproduksi dasar 𝑅 0 . Titik kesetimbangan bebas penyakit eksis tanpa syarat dan titik kesetimbangan endemi eksis jika 𝑅 0 > 1. Titik kesetimbangan bebas penyakit bersifat stabil  asimtotik  lokal  jika  nilai 𝑅 0 < 1 dan  titik  kesetimbangan  endemi  bersifat  stabil asimtotik lokal karena memenuhi kriteria Routh-Hurwitz.
 Pemberian  vaksinasi  dan screening memiliki  pengaruh  untuk  mengurangi  jumlah subpopulasi terinfeksi.
 Simulasi numerik yang dilakukan mendukung hasil analisis dinamik yang diperoleh.
5 Daftar Pustaka
[1] M. Nenad dan M. Olivera, What Every Woman Should Know about Cervical Cancer . Springer, 2001.
[2] F.X. Bosch, N. Munoz,dkk.,"Risk factors for cervical cancer in Colombia and Spain," Journal of cancer 52, pp. 750-758,1992.
[3] World  Health  Organization(WHO),  Human  Papillomavirus  (HPV)  and  Cervical  cancer. Diakses 12 Februari 2020, http://www.who.int/mediacentre/factsheets/fs380/en/ , 2016.
[4] T.  Malik  ,M.  Imran,  dan  R.  Jayaraman,  "Optimal  Control  with  Multiple  Human Papillomavirus," Journal of Theoretical Biology 393, pp. 179-193, 2016.
[5] O. Sharomi dan T. Malik, "A model to assess the effect of vaccine compliance on Human Papillomavirus  infection  and  cervical  cancer," Applied  Mathematical  Modelling 47,  pp. 528-550,2017.
[6] I.  Rasjidi,  "Epidemiologi  Kanker  Serviks," Indonesian  Journal  of  Cancer 3,  pp.  103- 108,2009.
[7] E.H.  Elbasha,  "Impact  of  Prophylactic  Vaccination  Againts  Human  Papillomavirus Infection," Contempory Mathematics 410, pp. 113-127,2006.
[8] L. Majed dan R. Lounes, "A SIS Model  for Human Papillomavirus Transmission," HAL archives-ouvertes 3, pp.15-23,2011.
[9] N. Van de Velde, M. Boily, M. D. Franco, E. Mayrand, dkk., "Population Level Impact of The Bivalent, Quadrivalent, and Nonavalent Human Papillomavirus Vaccines," Journal of Cancer 104, pp. 1712-1723, 2012.
[10]  F.  A.  Saldana,  I.  Korobeinikov,  dan  Barradas,  "Optimal  Control  Againts  The  Human Papillomavirus  Protection  Versus  Eradication  of  The  Infection," Abstract  and  Applied Analysis 10, pp. 1-13,2019.
[11]  Hasnawati, R. ratianingsih, dan J.W Puspita, "Analisis Kestabilan Model Matematika dan Penyebaran  Kanker  Serviks  Menggunakan  Kriteria  Routh-Hurwitz," Jurnal  Ilmiah Matematika dan Terapan 14, pp. 120-127, 2017.
[12]  N.  Puspitasari,  Y.  A.  Adi,  dan  R.S.  Winanda,  "Analisis  Kestabilan  Lokal  pada  Model Matematika  Kanker  Serviks  Akibat  Human  Papillomavirus," Jurnal  Ilmu  Alam  dan Teknologi Terapan , vol.1, no. 1, pp. 115-125, 2019.
[13]  F.  Brauer  dan  C.C.  Carlos,  " Mathematical  Models  in  Population  Biology  and Epidemiology," . Springer, 2010.
[14]  A. Panfilov, Qualitative Analysis of Differential Equations . Theoretical Biology,2004.
[15]  J. D. Murray, Mathematical Biology : An Introduction . Springer, 2002.