\documentclass{template-jurnal} %Bagian ini diedit oleh editor Jurnal Matematika UNAND
%\usepackage{array}
%\usepackage{booktabs}
\hyphenation{di-tulis-kan di-terang-kan}
%Perhatikan aturan penulisan dan ukuran huruf yang digunakan
\begin{document}

\markboth{Zulfatin Nafisah, Yudi Ari Adi } %Jika lebih dari dua penulis, tuliskan sebagai Nama Penulis Pertama dkk.
{Model SEIR dengan Pseudo-Recovery pada Kasus Tuberkulosis di Jawa Barat}

%%%%%%%%%%%%%%%%%%%%% Publisher's Area please ignore %%%%%%%%%%%%%%%
%
\catchline{}{}{}{}{}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\title{MODEL SEIR DENGAN PSEUDO-RECOVERY PADA KASUS TUBERKULOSIS DI JAWA BARAT}

\author{ZULFATIN NAFISAH, YUDI ARI ADI$^{}$\footnote{penulis korespondensi}\\}

\address{Program Studi Matematika, Fakultas Sains dan Teknologi Terapan, Universitas Ahmad Dahlan, Yogyakarta, Indonesia\\
email : \email{zulfatin2015015025@webmail.uad.ac.id , yudi.adi@math.uad.ac.id}}

\maketitle
\setcounter{page}{1} %bagian ini diedit oleh editor Jurnal Matematika UNAND

\begin{abstract}
\begin{center}
Diterima ..... \quad Direvisi ..... \quad Dipublikasikan ..... %tanggal-tanggal tersebut \textbf{dikosongkan} saja
\end{center}

\bigskip

\textbf{Abstrak}. %Dalam bahasa Indonesia
Tuberkulosis menjadi salah satu penyakit menular yang sampai saat ini masih sulit ditanggulangi. Di Indonesia, Provinsi Jawa Barat menjadi salah satu provinsi dengan kasus tertinggi yang memiliki tingkat prevalensi 63\%. Penelitian ini bertujuan untuk menganalisis model penyebaran penyakit tuberkulosis tipe SEIR (Susceptible-Exposed-Infected-Removed) dengan mempertimbangkan pengobatan dan adanya pemulihan semu (pseudo recovery). Model ini terdiri dari empat kelas yaitu, rentan, terpapar tetapi belum menularkan penyakit, terinfeksi dan dapat menularkan penyakit, dan sembuh semu. Data yang digunakan adalah data jumlah penderita penyakit tuberkulosis dari Dinas Kesehatan Provinsi Jawa Barat tahun 2011–-2022. Fungsi Lyapunov dan prinsip invarian LaSalle digunakan untuk menunjukkan bahwa titik keseimbangan stabil secara global, dan tuberkulosis akan bertahan jika angka reproduksi dasar lebih besar dari satu. Sebaliknya, penyakit akan hilang jika angka reproduksi dasar kurang dari satu.  Prosedur bifurkasi menggunakan teori manifold pusat digunakan untuk melakukan studi bifurkasi. Kondisi matematika memastikan terjadinya bifurkasi maju. Terakhir, simulasi numerik dilakukan untuk mendukung temuan teoretis..

\bigskip

\textbf{Abstract}. % Dalam bahasa Inggris
\textit{Tuberculosis is an infectious disease that is still difficult to overcome. In Indonesia, West Java Province is one of the provinces with the highest cases with a prevalence rate of 63\%. This study aims to analyse the SEIR (Susceptible-Exposed-Infected-Removed) model of tuberculosis disease spread by considering treatment and pseudo recovery. This model consists of four classes, namely, susceptible, exposed but not yet transmitting the disease, infected and can transmit the disease, and pseudo recovery. The data used is the number of tuberculosis patients from the West Java Provincial Health Office from 2011 to 2022. Lyapunov function and LaSalle invariance principle are used to show that the equilibrium point is globally stable, and tuberculosis will persist if the basic reproduction number is greater than one. Conversely, the disease will disappear if the basic reproduction number is less than one.  The bifurcation procedure using the central manifold theory is used to conduct bifurcation studies. Mathematical conditions ensure the occurrence of forward bifurcation. Finally, numerical simulations are performed to support the theoretical findings.}

\end{abstract}

\keywords{Analisis sensitifitas, Bifurkasi forward, Pseudo-recovery, Tuberkulosis}

\section{Pendahuluan}
Tuberkulosis merupakan salah satu penyakit menular yang disebabkan oleh bakteri Mycobacterium tuberculosis  \cite{1}. Data World Health Organization (WHO) menunjukkan bahwa sekitar 969 ribu orang di Indonesia terinfeksi tuberkulosis pada tahun 2021 dan sekitar 144 ribu orang meninggal karena tuberkulosis  \cite{2,3}. Tuberkulosis menempati peringkat ke-13 sebagai penyebab kematian terbanyak secara global dan merupakan penyakit menular kedua setelah COVID-19 dan HIV-AIDS  \cite{4}. Oleh karena itu, tuberkulosis menjadi masalah kesehatan utama di dunia.

Indonesia masih menjadi negara dengan tingkat kasus tuberkulosis tertinggi di dunia  \cite{6}. Jumlah kasus tuberkulosis yang tercatat berdasarkan diagnosa dokter dari 38 provinsi di Indonesia menunjukkan bahwa Provinsi Papua, Provinsi Banten, dan Provinsi Jawa Barat merupakan tiga daerah dengan prevalensi tertinggi, masing-masing mencapai 77\%, 76\%, dan 63\%  \cite{7}. Berdasarkan data tersebut Provinsi Jawa Barat menempati urutan ketiga dengan jumlah kasus tuberkulosis tertinggi di Indonesia. Menurut data Dinas Kesehatan Provinsi Jawa Barat, jumlah kasus tuberkulosis di provinsi tersebut meningkat dari tahun 2021 - 2023. Pada 2021, kasus baru mencapai 92.000 kasus, meningkat menjadi 159.000 pada 2022 \cite{8}. Penanganan tuberkulosis menjadi salah satu tujuan pembangunan berkelanjutan (SDGs)  \cite{9}.  Untuk memprediksi kemungkinan pencapaian tujuan tersebut, dapat digunakan model matematika.

Model matematika merupakan representasi formal yang dapat digunakan untuk menyederhanakan dan memahami masalah dalam kehidupan nyata. Salah satu contoh aplikasi model matematika adalah untuk mengetahui penyebaran penyakit menular di suatu daerah. Penyebaran penyakit menular dapat dimodelkan dengan menggunakan model deterministik atau stokastik. Model matematika tersebut diantaranya SIR, SEIR, SEIRS, SITR, SEIV. Setiap jenis model memiliki karakteristiknya masing-masing yang disesuaikan dengan jenis dan bentuk penyebaran penyakit menular yang diamati.

Pemodelan matematika terkait penyebaran tuberkulosis dapat dipelajari pada  \cite{10,11,12,13,14}. Syam, dkk  \cite{10} membahas tentang model SEIRS untuk penyebaran penyakit tuberkulosis di Kota Makassar. Mereka menggunakan data dari Kota Makassar dan menganalisis pola penyebaran penyakit tersebut dengan menggunakan model matematika. Chandra dan Roudhotillah  \cite{11} menganalisis stabilitas model penyebaran penyakit tuberkulosis dengan menggunakan pendekatan MSEITR. Mereka membahas implikasi hasil analisis stabilitas ini terhadap kontrol epidemiologi penyakit. Qomariyah dkk  \cite{12} melakukan analisis stabilitas pada model epidemiologi tuberkulosis dengan mempertimbangkan laju insidensi nonlinear dan efek pengobatan. Mereka menyoroti pentingnya efek pengobatan dalam strategi pengendalian penyakit. Selanjutnya Upadhyay dkk \cite{13} membahas tentang dinamika model epidemiologi SEIR dengan laju insidensi dan tingkat pengobatan nonlinier. Mereka menganalisis bagaimana faktor-faktor tersebut mempengaruhi penyebaran penyakit tuberkulosis. Sedangkan Qunia dkk  \cite{14} mengeksplorasi dampak tingkat kambuh pada model epidemiologi deterministik dengan pseudo-recovery. Mereka menyoroti pentingnya memperhitungkan tingkat kambuh dalam merancang kebijakan pengendalian penyakit. Artikel-artikel tersebut memberikan kontribusi penting dalam pemahaman tentang dinamika penyebaran penyakit tuberkulosis dan relevansinya dalam pengendalian penyakit tuberkulosis.

Pada artikel ini  dibahas model matematika tuberkulosis Susceptible-Exposed-Infected-Removed (SEIR) dengan Pseudo Recovery. Selanjutnya dari model yang terbentuk akan ditentukan analisis kestabilan lokal titik ekuilibrium bebas penyakit dan kestabilan global titik ekuilibrium bebas penyakit dan titik ekuilibrium endemik pada penyakit tuberkulosis. Dibahas juga analisis sensitifitas dan bifurkasi.


\section{Landasan Teori}
\noindent Pada bagian ini akan diberikan beberapa definisi dan teorema yang akan digunakan dalam pembahasan sistem persamaan diferensial, titik ekuilibrium, kriteria \textit{routh hurwitz}, dan bilangan reproduksi dasar.

\begin{definition} \cite{15}
 Diberikan sistem persamaan diferensial sebagai berikut :
\begin{equation}\label{eq1}
\frac{d}{dt}\textbf{x}=\dot{\textbf{x}}=\textbf{f(x)}
\end{equation}
dengan $\textbf{x}\in \mathbb{R}^{n}$ dan $\textbf{f}: E\rightarrow \mathbb{R}^{n}$ , E himpunan terbuka, $E\subseteq \mathbb{R}^{n}$.
 Fungsi $\textbf{x}(t)\in \mathbb{R}^{n}$ disebut solusi Sistem (\ref{eq1}) pada interval $I$  jika $\textbf{x}(t)$ diferensiabel pada $I$ dan untuk setiap $t\in I$ berlaku $\textbf{x}(t)\in E$ dan $\textbf{x}^{'}(t)=\textbf{f}(\textbf{x}(t))$, dengan $C(E)$ menyatakan himpunan semua fungsi kontinu pada $E$ dan $I$ interval terbuka.
\end{definition}

\begin{definition}\cite{15}
Titik $\overline{x}\ \in \ {\mathbb{R}}^n$\textit{ disebut titik ekuilibrium dari sistem persamaan diferensial }\eqref{eq1} \textit{jika memenuhi }$f\left(\overline{x}\right)=0.$\textit{}
\end{definition}

\begin{theorem} \cite{15}
Jika $f:\mathbb{R}^{n}\rightarrow \mathbb{R}^{n}$ memiliki turunan di $x_{0},$maka turunan parsial $\frac{\partial f_{i}}{\partial x_{j}},i,j=1,...,n$ ada pada $x_{0}$ dan untuk setiap $x\epsilon \mathbb{R}^{n}$ berlaku

\begin{equation} \label{eq2}
D\boldsymbol{f}(\boldsymbol{\bar{x}})\boldsymbol{x}=\sum_{i=1,j=1}^{n}\frac{\partial f_{i}}{\partial x_{j}}(\bar{\boldsymbol{x}})x_{j}
\end{equation}

matriks $D\boldsymbol{f}(\bar{\boldsymbol{x}})$ atau biasa ditulis dengan $J$
dinamakan matriks Jacobian.
\end{theorem}

Linearisasi sistem diberikan dalam definisi berikut.

\begin{definition} \cite{15}
Diberikan matriks Jacobian
$D\boldsymbol{f}(\bar{\boldsymbol{x}})$. Sistem linear $\dot{\boldsymbol{x}}=(D\boldsymbol{f}(\bar{\boldsymbol{x}}))\boldsymbol{x}$
disebut linearisasi dari sistem $\dot{\boldsymbol{x}}=\boldsymbol{f}(\boldsymbol{x})$
di sekitar titik $\bar{\boldsymbol{x}}$.
\end{definition}

\begin{definition} \cite{15}
Diberikan suatu sistem persamaan diferensial $\frac{dx}{dt}=f(x),$ dengan $x,f(x)\epsilon R$. Titik $\bar{x}\epsilon R$ disebut titik ekuilibrium jika dan hanya jika $f(\bar{\boldsymbol{x}})=0$.
\end{definition}

Kestabilan titik ekuilibrium Sistem \eqref{eq1} diberikan dalam definisi berikut.

\begin{definition} \cite{15}
Sistem linear dikatakan
\textbf{stabil} jika semua solusi pada sistem $t\rightarrow\infty$
dan \textbf{stabil asimtotik} jika semua solusi konvergen ke 0 dengan
$t\rightarrow\infty$.
\end{definition}


\begin{theorem} \cite{16}
 Misal persamaan karakteristik matriks $D\boldsymbol{f}(\bar{\boldsymbol{x}})$ pada sistem \eqref{eq2} dinyatakan sebagai

\begin{equation} \label{eq3}
\lambda^n+a_1{\lambda }^{n-1}+\dots +a_n=0
\end{equation}
dengan $a_n\neq 0.\ $Akar-akar persamaan karakteristik sistem $\left(2.1\right)$ memiliki bagian real negatif jika dan hanya jika
\[D_1=a_1>0,\ D_2=\ \left| \begin{array}{cc}
a_1 & a_3 \\
1 & a_2 \end{array}
\right|>0,\]


\noindent
\[D_3=\left| \begin{array}{ccc}
a_1 & a_3 & a_5 \\
1 & a_2 & a_4 \\
0 & a_1 & a_3 \end{array}
\right|>0,\dots ,\ D_k=\left| \begin{array}{c}
 \begin{array}{ccc}
a_1 & a_3 & a_5 \\
1 & a_2 & a_4 \\
0 & a_1 & a_3 \end{array}
 \begin{array}{ccc}
\dots  & \dots  & \dots  \\
\dots  & \dots  & \dots  \\
\dots  & \dots  & \dots  \end{array}
 \\
 \begin{array}{ccc}
0 & 1 & a_2 \\
\vdots  & \vdots  & \vdots  \\
0 & 0 & 0 \end{array}
\  \begin{array}{ccc}
\ \dots  & \dots  & \dots  \\
\dots  & \ddots  & \dots  \\
\dots  & \dots  & a_k \end{array}
 \end{array}
\right|>0,\]


\noindent dengan $k=1,2,\dots ,n$.

Untuk kasus tertentu, misalnya $n=4$, maka persamaan $\left(2.2\right)$  menjadi,

\begin{equation} \label{eq4}
{\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0.
\end{equation}

akar-akar persamaan \eqref{eq4} memiliki bagian real negatif jika dan hanya jika

\begin{enumerate}
\item  $a_1>0$

\item  $a_4>0$

\item  $a_1a_2-a_3>0$

\item  $\left(a_1a_2-a_3\right)a_3-a^2_1a_4>0$
\end{enumerate}
\end{theorem}

\begin{definition} \cite{17} \textit{Bilangan reproduksi dasar} ($R_0$) adalah jumlah rata-rata individu yang terinfeksi selama masa terinfeksinya dalam keseluruhan populasi rentan.
\end{definition}
Selanjutnya, persamaan kompartemen terinfeksi yang sudah dilinearkan dapat dituliskan sebagai berikut:
\begin{equation}
\dot{x} = (F - V)x
\end{equation}
Dengan $F$ dan $V$ merupakan matriks $n \times n$. $F_i$ adalah laju kemunculan infeksi baru yang menambah kelas terinfeksi, $V_i$ adalah laju perpindahan individu yang keluar dari kelas terinfeksi, $x_0$ adalah titik kesetimbangan bebas penyakit, dan $G$ merupakan matriks non negatif. Maka bilangan reproduksi dasar merupakan radius spektral dari matriks $G$ atau $\rho(G)$ yang dinyatakan dengan \cite{18,19,20}:
\begin{equation}
R_0 = \rho(G) = \rho(FV^{-1})
\end{equation}
dan bilangan reproduksi dasar diperoleh dari nilai eigen terbesar dari matriks $\left[ \frac{\partial F_i(x_0)}{\partial x_j} \right] \left[ \frac{\partial V_i(x_0)}{\partial x_j} \right]^{-1}$ dimana $F = \left[ \frac{\partial F_i(x_0)}{\partial x_j} \right]$ dan $V = \left[ \frac{\partial V_i(x_0)}{\partial x_j} \right]$.


\section{Pembahasan}
\subsection{Formulasi model}
Penyakit tuberkulosis merupakan salah satu penyakit menular yang dapat dimodelkan menggunakan model matematika epidemik. Untuk menyusun model epidemik tuberkuloasi ini, total populasi pada saat waktu $t$ disimbolkan dengan $N\left(t\right),\ $terdiri dari empat sub-populasi yaitu populasi \textit{susceptible} $\left(S\right)$ menyatakan jumlah individu yang rentan terhadap penyakit tuberkulosis, populasi \textit{exposed }$\left(E\right)$\textit{ }menyatakan individu yang terpapar tuberkulosis tetapi belum dapat menularkan penyakit, populasi \textit{infected} $\left(I\right)$ menyatakan jumlah individu yang terinfeksi dan dapat menularkan penyakit tuberkulosis, dan populasi \textit{removed} $\left(R\right)$ menyatakan jumlah individu yang telah sembuh. Formulasi model dalam penelitian ini menambahkan asumsi bahwa individu yang telah pulih terdapat kemungkinan kembali terinfeksi atau mengalami kekambuhan \textit{(pseudo recovery) }dan dapat menularkan penyakit. Diagram kompartemen model SEIR dengan \textit{pseudo recovery} penyakit tuberkulosis seperti pada Gambar 3.1.
\noindent
\begin{figure}[http]
\centering\includegraphics[width=9cm]{bagan}
\caption{Diagram kompartemen model SEIR dengan psudo-recovery} \label{gb1}
\end{figure}

Berdasarkan Gambar \ref{gb1}, diagram alur model SEIR dengan \textit{Pseudo Recovery} untuk dinamika penularan penyakit tuberkulosis; $\mathrm{\Lambda }$ adalah jumlah individu yang masuk karena adanya kelahiran dan migrasi, $\alpha SI\ $adalah kontak antara individu yang rentan dan terinfeksi, $\beta E$ adalah individu yang terinfeksi sepenuhnya dan dapat menularkan penyakit, $\frac{1}{1+\ \theta I}\ \ $adalah faktor penghambat \textit{(inhibitor)} terhadap laju kesembuhan. Faktor penghambat tersebut seperti biaya pengobatan  yang mahal, fasilitas kesehatan yang kurang memadai, persediaan obat yang terbatas, dan kebiasaan minum obat yang kurang teratur. Jumlah individu yang mati secara alami selama perjalanan penyakit dilambangkan dengan $\mu S;\mu E;\mu I;\mu R$, $\delta I$ adalah kematian yang disebabkan oleh penyakit di kelas \textit{infected} $\left(I\right)$, dan $\omega R$ adalah individu yang mengalami kekambuhan.

Berdasarkan Gambar \ref{gb1} diperoleh model dalam persamaan sebagai berikut :

\begin{equation} \label{eq31}
\begin{array} {rl}
 \frac{dS}{dt}= &  \Lambda-\left(\mu +\alpha I\right)S \\
\frac{dE}{dt}=  & \alpha SI-\left(\mu +\beta \right)E\\
 \frac{dI}{dt}= & \beta E-\left(\mu + \delta \right)I-\frac{\gamma I}{1+\ \theta I}+ \omega R \\
 \frac{dR}{dt}= & \frac{\gamma I}{1+\ \theta I}-\left(\mu +\omega \right)R
\end{array}
\end{equation}

yang memenuhi kondisi awal
\begin{equation} \label{eq32}
S\left(0\right),E\left(0\right),I\left(0\right),R\left(0\right)\ge 0
\end{equation}

dengan
\noindent %
\begin{tabular}{cl}
\emph{S} & : Individu rentan\tabularnewline
\emph{E} & : Individu yang terinfeksi tuberkulosis tetapi tidak dapat menularkan penyakit\tabularnewline
\emph{I} & : Individu yang terinfeksi dan dapat menularkan penyakit tuberkulosis\tabularnewline
\emph{R} & : Individu yang telah sembuh\tabularnewline
$\Lambda$ &: Jumlah individu yang masuk (kelahiran, migrasi)\tabularnewline
$\alpha $ &: Laju individu rentan menjadi terinfeksi\tabularnewline
$\beta $ &: Laju penularkan penyakit tuberkulosis\tabularnewline
$\gamma $ &: Laju kesembuhan dari infeksi yang aktif\tabularnewline
$\theta $ &: Faktor \textit{inhibitor }(penghambat) terhadap penularan penyakit tuberkulosis\tabularnewline
$\mu $ &: Laju kematian alami\tabularnewline
$\delta $ &: Laju kematian karena terinfeksi tuberkulosis pada setiap individu\tabularnewline
$\omega $ &: Laju kekambuhan penyakit tuberkulosis\tabularnewline
\end{tabular}

\noindent
Selanjutnya untuk menunjukkan bahwa variabel dalam model bernilai positif, maka akan dibuktikan bahwa sistem \eqref{eq31} mempunyai solusi positif dan terbatas dalam Teorema 3.1.

\begin{theorem} Himpunan $\mathrm{\Omega }=\left\{\left(S,E,I,R\right)\in \ {\mathbb{R}}^4_+\ :S+E+I+R\ \le \ \frac{\mathrm{\Lambda }}{\mu }\right\}$ merupakan invarian positif untuk sistem \eqref{eq31} dengan kondisi awal \eqref{eq32}.
\end{theorem}

\noindent \textbf{Bukti.}
Berdasarkan persamaan \eqref{eq31} diperoleh
\[\frac{dN}{dt}=\frac{dS}{dt}+\frac{dE}{dt}+\frac{dI}{dr}+\frac{dR}{dt}\]
\[=\ \mathrm{\Lambda }-\left(\mu +\alpha I\right)S+\ \alpha SI-\left(\mu +\beta \right)E+\ \beta E-\left(\mu +\ \delta \right)I-\frac{\gamma I}{1+\ \theta I}+\ \omega R\ +\ \frac{\gamma I}{1+\ \theta I}-\left(\mu +\omega \right)R\]
\[=\ \ \mathrm{\Lambda }\mathrm{-}\mu S-\ \mu E-\ \mu I-\ \mu R-\ \delta I\mathrm{\ }\]
\[=\mathrm{\Lambda }\mathrm{-}\mu \left(S+E+I+R\right)-\delta I\]
\[=\ \mathrm{\Lambda }\mathrm{-}\mu N-\ \delta I\]
sehingga diperoleh
\[\frac{dN}{dt}=\ \mathrm{\Lambda }\mathrm{-}\mu N-\ \delta I\ \le \ \mathrm{\Lambda }\mathrm{-}\mu N.\]

Selanjutnya dengan mengintegralkan pertidaksamaan di atas diperoleh penyelesaian menggunakan kondisi awal $t=0,\ N_0\ $sebagai berikut :
\[N\left(t\right)\le N_0e^{-\mu t}+\frac{\mathrm{\Lambda }}{\mu }\]
Dengan demikian, jika kondisi awal berada dalam himpunan $\mathrm{\Omega }$, maka solusi sistem akan selalu berada dalam $\mathrm{\Omega }$ untuk $t\to \ \infty .$ Sebaliknya, jika kondisi awal berada di luar himpunan $\mathrm{\Omega }$, maka solusi sistem akan cenderung mendekati himpunan $\mathrm{\Omega }$ untuk $t\to \ \infty .$ Oleh karena itu, himpunan $\mathrm{\Omega }$ merupakan invarian positif.

\subsection{Titik ekuilibrium dan analisis kestabilan}
Titik ekuilibrium Sistem \eqref{eq31} terjadi pada saat $\left(\frac{dS}{dt},\frac{dE}{dt},\frac{dI}{dt},\frac{dR}{dt}\right)=\left(0,0,0,0\right)$ Diperoleh titik ekuilibrium bebas penyakit
\[E_0=\left(S^0,E^0,I^0,R^0\right)=\left(\frac{\mathrm{\Lambda }}{\mu }\ ,0,0,0\right).\]

Berikutnya ditentukan bilangan reproduksi dasar. Bilangan reproduksi dasar merupakan bilangan yang menyatakan nilai rata-rata dari individu pada populasi rentan yang terinfeksi secara langsung oleh individu pada populasi terinfeksi. Metode \textit{next generation matrix }digunakan untuk memperoleh bilangan reproduksi dasar, yaitu $G=FV^{-1}$. Dalam model epidemiologi ini, individu yang mengalami kekambuhan penyakit dikelompokkan ke matriks $\mathcal{F}$. Matriks $\mathcal{F}$ merupakan kompartemen $E$, $I,\ $dan $R$ yang merupakan kelas individu yang telah pulih dari infeksi, tetapi masih memiliki kemungkinan kambuh kembali. Sedangkan matriks $\mathcal{V}$ merupakan transisi yang tersisa yaitu laju kelahiran, laju kematian, penyebaran penyakit, dan laju pemulihan. $R_0$ didefinisikan sebagai radius \textit{spektral} dari \textit{next generation matrix }$FV^{-1}$\textit{ }dengan $F$ dan $V$ adalah matriks jacobian dari $\mathcal{F}$ dan $\mathcal{V}$ pada titik ekuilibrium bebas penyakit. Oleh karena itu dari persamaan \eqref{eq31}diperoleh,
\[\mathcal{F}\left[E,I,R\right]=\ \left[ \begin{array}{c}
 \begin{array}{c}
{\mathcal{F}}_1 \\
{\mathcal{F}}_2 \end{array}
 \\
{\mathcal{F}}_3 \end{array}
\right]=\left[ \begin{array}{c}
 \begin{array}{c}
\alpha SI \\
0 \end{array}
 \\
0 \end{array}
\right]\]
\textbf{}
\[\mathcal{V}\left[E,I,R\right]=\left[ \begin{array}{c}
{\mathcal{V}}_1 \\
 \begin{array}{c}
{\mathcal{V}}_2 \\
{\mathcal{V}}_3 \end{array}
 \end{array}
\right]=\left[ \begin{array}{c}
\left(\mu +\beta \right)E \\
 \begin{array}{c}
-\beta E+\left(\mu +\delta \right)I+\frac{\gamma I}{1+\ \theta I}-\ \omega R \\
-\frac{\gamma I}{1+\ \theta I}+\left(\mu +\omega \right)R \end{array}
 \end{array}
\right]\]

maka diperoleh bilangan reproduksi dasar ${(R}_0),\ $yaitu
\[R_0= \frac{\alpha \Lambda \beta \left(\mu +\omega \right)}{\mu \left(\mu + \beta \right)\left(\mu \left(\mu + \delta +\gamma \right)+\omega \left(\mu +\delta \right)\right)}\]

Kestabilan titik ekuilibrium bebas penyakit diberikan dalam teorema berikut.
\begin{theorem} \label{th32}
Titik ekuilibrium bebas penyakit $\left(E_0\right)$ stabil asimtotik lokal jika $R_0<1.$
\end{theorem}

\noindent \textbf{Bukti. } Matriks Jacobian sistem \eqref{eq31} di titik ekuilibrium bebas penyaki $E_0=\left(\frac{\mathrm{\Lambda }}{\mu },0,0,0\right)$ adalah
\[
J(S^0, E^0, I^0, R^0) = \begin{bmatrix}
    -\mu & 0 & -\alpha \left( \frac{\Lambda}{\mu} \right) & 0 \\
    0 & -( \mu + \beta) & \alpha \left( \frac{\Lambda}{\mu} \right) & 0 \\
    0 & \beta & -(\mu + \delta + \gamma) & \omega \\
    0 & 0 & \gamma & -(\mu + \omega)
\end{bmatrix} \]

dan persamaan karakteristiknya adalah
\[(\lambda +\mu) \left[\lambda +\left(\mu +\beta \right)\left[\left(\lambda +\mu +\delta +\gamma \right)\left(\lambda +\mu +\omega \right)-\omega \gamma \right]+\beta \left(\frac{-\alpha \mathrm{\Lambda }}{\mu }\left(\lambda +\mu +\omega \right)\right)\right]=0\ \]

sehingga diperoleh nilai eigen
\[{\lambda }_1=\ -\mu <0\]

dan tiga nilai eigen lainnya merupakan akar persamaan kubik
\begin{equation} \label{eq33}
{\lambda }^3+a_2{\lambda }^2+a_1\lambda +a_0=0,
\end{equation}
dengan,
\begin{eqnarray*}
a_2 &=&\left(\mu +\omega \right)+\left(\mu +\delta +\gamma \right)+\left(\mu +\beta \right)>0 \nonumber\\
  a_1 &=&\ \left(\mu +\beta \right)\left(\ \mu +\delta +\gamma \right)+\left(\mu +\beta \right)\left(\mu +\omega \right)+\left(\ \mu +\delta +\gamma \right)\left(\delta +\gamma \right)-\left(\omega \gamma -\frac{\beta \alpha \mathrm{\Lambda }}{\mu }\right) \\
a_0&=&\left(\mu +\beta \right)\left(\mu \left(\mu +\delta +\gamma \right)+\omega \left(\mu +\delta \right)\right)\left(1-R_0\right). \nonumber\\
\end{eqnarray*}
Menurut kriteria \textit{Routh-Hurwitz } persamaan \eqref{eq33} memiliki akar-akar dengan bagian real negatif jika $a_2>0,a_0>0,\ $ dan $a_2a_1-a_0>0.$

\begin{eqnarray*}
a_2 a_1 - a_0 &=& \left[ (\mu + \omega) + (\mu + \delta + \gamma) + (\mu + \beta) \right]
\left[ (\mu + \beta)(\mu + \delta + \gamma) + (\mu + \beta)(\mu + \omega)\right] \\ \nonumber
&+&(\mu + \delta + \gamma)(\delta + \gamma) - (\omega \gamma - \frac{\beta \alpha \Lambda}{\mu})
- (\mu + \beta) (\mu + \delta + \gamma)(\mu + \omega) - \omega \gamma \\ \nonumber
&-& \frac{\beta \alpha \Lambda}{\mu} (\mu + \omega)  \\ \nonumber
&=& \left[ (\mu + \omega) + (\mu + \delta + \gamma) + (\mu + \beta) \right]
 [(\mu + \beta) ( (\mu + \delta + \gamma) + (\mu + \omega) ) \\ \nonumber
 &+& (\mu + \delta + \gamma)(\delta + \gamma) + \frac{\beta \alpha \Lambda}{\mu} - \omega \gamma
-  (\mu + \beta) (\mu (\mu + \delta + \gamma) + \omega (\mu + \delta) \\ \nonumber
&-& \frac{\beta \alpha \Lambda}{\mu} (\mu + \omega) ) \\ \nonumber
&=& \left[ (\mu + \beta) \left[ (\mu + \delta + \gamma) + (\mu + \omega) \right] + (\mu + \delta + \gamma)(\delta + \gamma) + \frac{\beta \alpha \Lambda}{\mu} - \omega \gamma \right] A \\ \nonumber
&-& B,
\end{eqnarray*}
dengan
\begin{eqnarray*}
A &=& (\mu + \omega) + (\mu + \delta + \gamma) + (\mu + \beta) \\ \nonumber
B &=& (\mu + \beta)\left[ \mu (\mu + \delta + \gamma) + \omega (\mu + \delta) \right] \left(1 - \frac{\beta \alpha \Lambda (\mu + \omega)}{(\mu + \beta) [\mu (\mu + \delta + \gamma) + \omega (\mu + \delta)]} \right)
\end{eqnarray*}
Untuk menunjukkan semua nilai eigen memiliki bilangan real negatif digunakan kriteria \textit{Routh Hurwitz} yaitu $a_2>0,\ \ a_0>0{,\ \ a}_2a_1-a_0>0$, maka dapat disimpulkan bahwa titik ekuilibrium bebas penyakit ($E_0)$ stabil asimtotik lokal dengan syarat $R_0<1.$
Selanjutnya dapat diperoleh titik ekuilibrum endemik yag dimyatakan dalam teorema berikut.

\begin{theorem} \label{th33}
Sistem \eqref{eq31} mempunyai:
\begin{enumerate}
\item  satu titik ekulibrium endemik jika $R_0 > 1$
\item tidak ada titik ekuilibrum endemik jika $R_0<1$.
\end{enumerate}
\end{theorem}
\noindent \textbf{Bukti.}
Titik ekuilibrium tercapai pada saat $\left(I\neq 0\right)$. Dari sistem \eqref{eq1} diperoleh

\begin{equation} \label{eq34}
\begin{array} {rl}
\Lambda-\left(\mu +\alpha I^*\right)S^*=0\\
\alpha S^*I^*-\left(\mu +\beta \right)E^*=0  \\
 \beta E^*-\left(\mu +\delta \right)I^*-\ \frac{\gamma I^*}{1+\ \theta I^*}+\ \omega R^*=0\\
\frac{\gamma I^*}{1+\ \theta I^*}-\left(\mu +\omega \right)R^*=0
\end{array}
\end{equation}

sehingga $S^*=\frac{\Lambda}{\mu+\alpha I^*}$, $E^*=\ \frac{\alpha \mathrm{\Lambda }I^*}{\left(\mu +\beta \right)\left(\mu +\alpha I^*\right)}$, dan $R^*=\ \frac{\gamma I^*}{\left(1+\ \theta I^*\right)\left(\mu +\omega \right)}$. Selanjutnya substitusi ke  ke dalam persamaan ketiga pada \eqref{eq34} diperoleh,
\[a_2I^2+a_1I+a_0=0\]
dengan
\begin{equation} \label{eq34}
\begin{array} {rl}
{\ a}_2= & -\left(\mu +\omega \right)\left(\mu +\delta \right)\left(\mu +\beta \right)\alpha \theta <0 \\
a_1 = &\ \left(\mu +\omega \right)\beta \alpha \mathrm{\Lambda }\theta -\left(\mu +\beta \right)\left(\left(\mu +\omega \right)\left(\mu +\beta \right)\left(\mu +\delta +\gamma \right)+\mu \alpha \gamma \right)\\
a_0=& \left(R_0-1\right)\mu \left(\mu +\beta \right)\left(\mu \left(\mu +\delta +\gamma \right)+\omega \left(\mu \mathrm{+}\mathrm{\delta }\right)\right)<0,\ $jika $R_0<1 \\
 = & \left(R_0-1\right)\mu \left(\mu +\beta \right)\left(\mu \left(\mu +\delta +\gamma \right)+\omega \left(\mu \mathrm{+}\mathrm{\delta }\right)\right)>0,\ $jika $R_0>1.
\end{array}
\end{equation}

\noindent Berdasarkan aturan tanda Descartes, maka diperoleh jumlah perubahan tanda pada titik ekuilibrium endemik $\left(I^*\right)$ sebagai berikut :
\noindent
\begin{table}[htbp]
\begin{center}
\begin{small}
\caption{Titik ekuilibrium endemik}\label{tab1}
\begin{tabular}{|c|c|c|c|c|}
\hline	 & $a_2$ & $a_1$	& $a_0$ & Banyaknya akar positif	\\
\hline $R_0<1$&	 $<0$ &	$<0$	&	$<0$	& $0$ 	\\
\hline $R_0>1$&	 $<0$ &	$<0$	&	$>0$	& $1$	\\
\hline								
\end{tabular}
\end{small}

\end{center}
\end{table}

Dari Tabel 1 terlihat bahwa jika $R_0<1$ maka tidak ada titik ekuilibrium endemik, sedangkan jika $R_0>1$ hanya terdapat satu titik ekuilibrium endemik. Bukti selesai.

Kestabilan titik ekuilibrium endemik dinyatakan dalam teorema berikut.
\begin{theorem} \label{th34}
Titik ekuilibrium endemik \(E_1 = (S^*, E^*, I^*, R^*)\) pada sistem \eqref{eq31} stabil asimtotik global jika \(R_0 > 1\).
\end{theorem}
\section*{Bukti}

Berikut diberikan Fungsi Lyapunov tipe Volterra:
\begin{align*}
L = & \left( S - S^* - S^* \ln \frac{S}{S^*} \right) + \left( E - E^* - E^* \ln \frac{E}{E^*} \right) \\
    & + \frac{\mu + \beta}{\beta} \left( I - I^* - I^* \ln \frac{I}{I^*} \right) \\
    & + \frac{\omega (\mu + \beta)}{\beta (\mu + \omega)} \left( R - R^* - R^* \ln \frac{R}{R^*} \right)
\end{align*}

Turunan Fungsi Lyapunov diperoleh sebagai berikut:
\begin{align*}
\frac{dL}{dt} = & S - \frac{S S^*}{S} + E - \frac{E E^*}{E} + \frac{\mu + \beta}{\beta} \left( I - \frac{I I^*}{I} \right) \\
                & + \frac{\omega (\mu + \beta)}{\beta (\mu + \omega)} \left( R - \frac{R R^*}{R} \right) \\
              = & S \left( 1 - \frac{S^*}{S} \right) + E \left( 1 - \frac{E^*}{E} \right) + \frac{(\mu + \beta)}{\beta} I \left( 1 - \frac{I^*}{I} \right) \\
                & + \frac{\omega (\mu + \beta)}{\beta (\mu + \omega)} R \left( 1 - \frac{R^*}{R} \right) \\
              = & \Lambda \left( 1 - \frac{S^*}{S} \right) - \mu S \left( 1 - \frac{S^*}{S} \right) + \alpha I S^* - \alpha I S \frac{E^*}{E} \\
                & + (\mu + \beta) E^* - \frac{(\mu + \beta) (\mu + \delta) I}{\beta} - \frac{(\mu + \beta) \gamma I}{\beta (1 + \theta I)} \\
                & - (\mu + \beta) E \frac{I^*}{I} + \frac{(\mu + \beta) (\mu + \delta) I^*}{\beta} + \frac{(\mu + \beta) \gamma I^*}{\beta (1 + \theta I)} \\
                & - \frac{(\mu + \beta) \omega R I^*}{\beta I} + \frac{\omega (\mu + \beta)}{\beta (\mu + \omega)} \left( \frac{\gamma I}{1 + \theta I} \right) \\
                & - \frac{\omega (\mu + \beta)}{\beta (\mu + \omega)} \left( \frac{\gamma I}{1 + \theta I} \right) \frac{R^*}{R} + \frac{\omega (\mu + \beta)}{\beta} R^*
\end{align*}
Dari sistem (3.1) diperoleh
\begin{align*}
\Lambda &= (\mu + \alpha I^*) S^*, \\
\mu + \beta &= \frac{\alpha S^* I^*}{E^*}, \\
\mu + \delta &= \frac{\beta E^* - \frac{\gamma I^*}{1 + \theta I^*} + \omega R^*}{I^*}, \\
\mu + \omega &= \frac{\gamma I^*}{(1 + \theta I^*) R^*}.
\end{align*}

sehingga
\begin{align*}
\frac{dL}{dt} &= \mu S^* \left[ 2 - \frac{S^*}{S} - \frac{S}{S^*} \right] + \alpha S^* I^* \left[ 3 - \frac{S^*}{S} - \frac{SI}{S^* I^*} \frac{E^*}{E} - \frac{EI^*}{E^* I} \right] \\
&\quad + \frac{\gamma I^*}{E^* \beta (1 + \theta I)} \left[ 1 - \frac{I}{I^*} \right] + \frac{\omega R^*}{E^* \beta} \left[ 1 - \frac{R}{R^*} \right] \\
&\quad + \frac{\alpha S^* I^* \omega R^*}{E^* \beta} \left[ 1 - \frac{RI^*}{R^* I} \right] + \frac{\omega \alpha S^* (1 + \theta I^*) R^* I}{E^* \beta (1 + \theta I)} \left[ 1 - \frac{R^*}{R} \right]
\end{align*}

Dengan menggunakan AM-GM (Arithmetic Mean – Geometric Mean)\cite{21}, diperoleh:
\begin{align*}
2 - \frac{S^*}{S} - \frac{S}{S^*} &\leq 0, \\
3 - \frac{S^*}{S} - \frac{SI}{S^* I^*} \frac{E^*}{E} - \frac{EI^*}{E^* I} &\leq 0, \\
1 - \frac{I}{I^*} & \leq 0, \\
1 - \frac{R}{R^*} & \leq 0,\\
1 - \frac{R I^*}{R^* I} & \leq 0,\\
1 - \frac{R^*}{R} &  \leq 0.\\
\end{align*}

Sehingga $
\frac{dL}{dt} \leq 0 \text{ dengan } \frac{dL}{dt} = 0 \text{ jika } S = S^*, E = E^*, I = I^*, R = R^*$. Dengan demikian diperoleh $
\frac{dL}{dt} \leq 0$. Berdasarkan teorema Lyapunov-LaSalle’s invariance, invarian subset terbesar dimana \(\frac{dL}{dt} = 0\) yaitu \(E_1 = (S^*, E^*, I^*, R^*)\) dan dapat disimpulkan bahwa titik ekuilibrium endemik \(E_1\) stabil asimtotik global jika \(R_0 > 1\). Bukti selesai.

Dari teorema ini dapat disimpulkan bahwa penyakit tuberkulosis akan menyebar kapan pun ketika \(R_0 > 1\) berapapun jumlah awal individu yang menularkan penyakit tersebut dalam populasi.

\section{Simulasi numerik}
\subsection{Estimasi parameter}
Pada bagian ini model akan diselesaikan secara numerik. Sebagai analisis awal untuk menentukan parameter model epidemi tuberkulosis, beberapa parameter diasumsikan berdasarkan data tuberkulosis yang diperoleh dari website Dinas Kesehatan Provinsi Jawa Barat pada tahun $2011-2022$. Usia harapan hidup rata-rata penduduk Jawa Barat dari tahun $2011-2022$ adalah 66 tahun, maka \(\mu = \frac{1}{66} = 0.01515\). Badan Pusat Statistik (BPS) menyatakan jumlah penduduk Jawa Barat di tahun 2011 yaitu 89.174.917 jiwa sehingga untuk menghitung laju recruitment (\(\Lambda\)) yaitu, tingkat kematian alami \(\times\) total subpopulasi, yang pada kasus ini menghasilkan \(\Lambda = 1.351.000\) \cite{22}. Capaian keberhasilan pengobatan di Provinsi Jawa Barat pada tahun 2020 mencapai 83\% dari target yang ditetapkan sebesar 90\%. Angka tersebut turun menjadi 74\% dari tahun 2021 hingga triwulan III tahun 2022. Berdasarkan hal tersebut, target penurunan insidensi tuberkulosis untuk tahun 2024 telah ditetapkan sebesar 190 per 100.000 penduduk \cite{23}. Untuk parameter lainnya, menggunakan metode least square melalui minimalisasi fungsi objektif kuadratik.

\[
\text{SSE} = \sum_{j=1}^n (I_{a_j} - I_{m_j})^2
\]

Dimana \(I_{a_j}\), \(I_{m_j}\), \(n\) mewakili masing-masing jumlah kasus infeksi tuberkulosis aktual, jumlah kasus tuberkulosis yang diprediksi, dan jumlah total kasus. Untuk keperluan fitting data, penulis menghitung total kasus tuberkulosis yang dikonfirmasi sebagai jumlah total individu yang terpapar (\(E\)) dan individu yang terinfeksi (\(I\)). Nilai parameter yang paling sesuai dengan kasus tuberkulosis diberikan pada Tabel \ref{tab2}. Hasil fitting antara model dengan data penderita tuberkulosis di Jawa Barat dari tahun $2011-2022$ ditunjukkan pada Gambar \ref{gb2}. Tingkat akurasi berdasarkan nilai \textit{mean absolute percentage error} (MAPE) sebesar 8.5974\% yang menunjukkan bahwa hasil estimasi parameter yang diperoleh akurat.

\begin{table}[htbp]
\begin{center}
\begin{small}
\caption{Nilai parameter}\label{tab2}
\begin{tabular}{|c|c|c|}
\hline	Parameter & nilai	& sumber 	\\
\hline \(\Lambda\)  & 1.351.000 & \cite{22} \\
\hline \(\mu\) &  0.01515 & \cite{22} \\
\hline\(\alpha\)   & \(1.0563 \times 10^{-8}\) & Fitting \\
\hline\(\beta\)  & 0.0010921 & Fitting \\
\hline\(\delta\)  & 0.0072373 & Fitting \\
\hline\(\gamma\)  & 0.0083693 & Fitting \\
\hline\(\theta\)  & 0.0012755 & Fitting \\
\hline\(\omega\)  & 0.031165 & Fitting \\
\hline
\end{tabular}
\end{small}
\end{center}
\end{table}


\begin{figure}[http]
\centering\includegraphics[width=9cm]{ep1}
\caption{Fitting model dengan data tuberkulosis di Jawa Barat} \label{gb2}
\end{figure}

Selanjutnya dengan menggunakan nilai-nilai parameter pada Tabel \ref{tab2} diperoleh $R_0=2.5208 >1$, yang menunjukkan bahwa titik ekuilibrium endemik stabil global. Gambar \ref{gb3} menunjukkan subpopulasi terpapar (E),  subpopulasi terinfeksi (I), dan subpopulasi pulih (R) menggunakan empat nilai awal yang berbeda akan menuju ke titik ekuilibrium  $E_1=(S^*,E^*,I^*,R^* )=(3.1521 ,5.3776,2.6232,141.63)$. Hal ini menunjukkan bahwa titik ekuilibrium endemik $E_1$ bersifat stabil asimtotik. Dengan demikian ketika bilangan reproduksi dasar $R_0>1$, maka dengan kondisi awal apapun penyakit akan tetap ada dalam populasi dan penyebarannya terus terjadi.
\begin{figure}[http]
\centering\includegraphics[width=11cm]{end1}
\caption{Trayektori fase menggunakan nilai parameter pada Tabel \ref{tab2}} \label{gb3}
\end{figure}

 Berdasarkan Teorema \ref{th33} dan Teorema \ref{th34} dapat disimpulkan terjadi bifurkasi forward ketika $R_0=1$. Hal ini dinyatkan dalam teorema berikut.

\begin{theorem} \label{th41}
 Sistem (\ref{eq31})mengalami bifurkasi forward di  $R_0 =1$, yang bersesuaian dengan $\alpha=\alpha^*=\frac{\mu (\beta+\mu)(\mu(\gamma+\delta+\mu)+\omega(\mu+\delta))}{\Lambda \beta (\mu+\omega)}$.
\end{theorem}
\noindent
\textbf{Bukti.} Matriks Jacobian sistem (\ref{eq31}) di $E_0$ dengan $\alpha=\alpha^*$ adalah
\[J(E_0,\alpha^*)=\left[\begin{array}{cccc}
-\mu& 0 & -\alpha^*\frac{\Lambda}{\mu} & 0\\
0 &-(\beta+\mu)&\alpha^*\frac{\Lambda}{\mu}&0\\
0& \beta &-(\gamma+\theta+\mu) &\omega\\
0 & 0 & \gamma & -(\mu+\omega)
\end{array}\right],\]
mempunyai tiga nilai eigen real negatif dan satu nilai eigen nol. Dengan demikian, $E_0$ adalah titik ekuilibrium non-hyperbolic, sehingga sistem (\ref{eq31}) mengalami bifurkasi saat $R_0=1$, yaitu di $\alpha=\alpha^*$. Berikutnya akan ditentukan vektor eigen yang bersesuaian dengan nilai eigen nol.

Eigenvektor kanan, $\textbf{w}$ dan eigenvektor kiri, $\textbf{v}$ yang bersesuaian dengan nilai eigen nol diperoleh dari  $J(E_0,\alpha^*)\textbf{w}=0$ dan $\textbf{v}J(E_0,\alpha^*)=0$.
vektor eigen kanan
\[\textbf{w}=\left(\frac{-\alpha^*\Lambda}{\mu^2}w_3, \frac{\alpha^*\Lambda}{\mu(\beta+\mu)}w_3,w_3,\frac{\gamma}{\mu+\omega}w_3\right)^T\]
dan
\[\textbf{v}=\left(0,\frac{\beta}{\beta+\mu}v_3,v_3,\frac{\omega}{\mu+\omega}v_3\right).\]

Agar $\textbf{v.w}=1$, maka $w_3$ dan $v_3$ adalah
\[w_3=\frac{1}{\alpha^*\beta\Lambda(\omega+\mu)^2+\mu(\mu+\beta)^2(\mu+\omega)^2+\omega \gamma \mu(\mu+\beta)^2}\]
dan
\[v_3=\mu(\beta+\mu)^2(\mu+\omega)^2.\]

 Selanjutnya, arah bifurkasi dapat dilihat dari koefisien bifurkasi $a$ dan $b$, dengan
 \[a=\sum\limits_{k,i,j=1}^{4}{{{v}_{k}}{{w}_{i}}{{w}_{j}}\frac{{{\partial }^{2}}{{f}_{k}}}{\partial {{x}_{i}}\partial {{x}_{j}}}\left( {{E}_{0}},\alpha^* \right),\qquad } b=\sum\limits_{k,i=1}^{4}{{{v}_{k}}{{w}_{i}}\frac{{{\partial }^{2}}{{f}_{k}}}{\partial {{x}_{i}}\partial \alpha}\left( {{E}_{0}},\alpha^* \right)}.\]

Dengan memperhatikan komponen dari eigen vektor kanan dan eigen vektor kiri $\textbf{w}$ dan $\textbf{v}$ dan turunan parsial kedua dari $f_i, i=1,2,3,4$ dari sistem (\ref{eq31}) diperoleh
\begin{eqnarray*}
 a&=&v_2w_1w_3\frac{\partial^2f_2}{\partial {S}\partial {I}}(E_0,\alpha^*) +{{v}_{2}}{{w}_{3}}{{w}_{1}}\frac{{{\partial }^{2}}{{f}_2}}{\partial {I}\partial {S}}({{E}_{0}},\alpha^*)\\
    &=&-2\frac{\alpha^{*3} \Lambda^2 v_2 w_3^{2}}{\mu^3}\\
  \end{eqnarray*}
  dan
 \begin{eqnarray*}
 b&=&v_2w_3\frac{\partial^2f_2}{\partial {I}\partial {\alpha}}(E_0,\alpha^*) \\
  &=&\frac{\alpha^* \Lambda \beta v_3 w_3}{\mu(\beta+\mu)}.
  \end{eqnarray*}

Mengingat $v_3, w_3$, dan semua nilai parameter sistem (\ref{eq31}) adalah positif, maka $a<0$ dan $b>0$ sehingga menurut Teeorem 4 di dalam \cite{24}, disimpulkan bahwa sistem (\ref{eq31}) mempunyai bifurkasi forward di $R_0=1$, atau ekuivalen dengan  $\alpha=\alpha^*=\frac{\mu (\beta+\mu)(\mu(\gamma+\delta+\mu)+\omega(\mu+\delta))}{\Lambda \beta (\mu+\omega)}$. Bukti selesai.

\subsection{Analisis sensitivitas}
Analisis sensitivitas digunakan untuk menentukan parameter yang paling berpengaruh terhadap bilangan reproduksi dasar $R_0$ yang menunjukkan terjadinya penyebaran penyakit tuberkulosis di dalam populasi. Oleh karena itu, analisis sensitivitas penting dilakukan untuk mengendalikan penyebaran penyakit tuberkulosis. Indeks sensitivitas merupakan rasio perubahan relatif dalam $R_0$ terhadap perubahan relatif dalam parameter $p$ seperti berikut \cite{25}.
\[
C_{R_0}^p = \frac{\partial R_0}{\partial p} \times \frac{p}{R_0}.
\]

Indeks sensitivitas parameter terhadap $R_0$  disajikan dalam Tabel \ref{tab3}.
\begin{table}[htbp]
\begin{center}
\begin{small}
\caption{Nilai parameter}\label{tab3}
\begin{tabular}{|c|c|}
\hline	Parameter & Indeks sensitivitas 	\\
\hline \(\mu\)  &$-2.6090$ \\
\hline \(\delta\) &  $-0.2880$ \\
\hline\(\gamma\)   & $-0.1089$ \\
\hline\(\beta\)  & 0.9327 \\
\hline\(\omega\)  & 0.0733 \\
\hline\(\alpha\)  & 1 \\
\hline
\end{tabular}
\end{small}
\end{center}
\end{table}

Indeks sensitivitas parameter terhadap $R_0$  disajikan dalam Tabel \ref{tab3}. Berdasarkan Tabel \ref{tab3} diperoleh indeks sensitivitas parameter yang bernilai positif adalah  \(\alpha\), \(\beta\), dan \(\omega\). Hal ini menunjukkan bahwa jika salah satu dari parameter \(\alpha\), \(\beta\), dan \(\omega\) dinaikkan sementara parameter yang lain konstan, maka nilai \(R_0\) akan meningkat sehingga menyebabkan penyebaran penyakit tuberkulosis akan selalu terjadi. Parameter yang memiliki nilai indeks sensitivitas negatif adalah \(\mu\), \(\delta\), dan \(\gamma\). Hal ini menunjukkan bahwa jika salah satu dari parameter \(\mu\), \(\delta\), dan \(\gamma\) dinaikkan sementara parameter yang lain konstan, maka nilai \(R_0\) akan menurun sehingga menyebabkan penyebaran penyakit tuberkulosis juga menurun.

Indeks sensitivitas \(\mu\) memiliki tanda negatif yang berarti semakin besar tingkat kematian semakin berkurang \(R_0\). Namun, tidak mungkin untuk meningkatkan tingkat kematian alami meskipun memiliki indeks sensitivitas yang paling luas. Parameter lain yang memiliki indeks sensitivitas besar adalah  \(\alpha\) dan \(\beta\) yang memiliki indeks sensitivitas 1 dan 0.9327. Oleh karena itu, untuk mengurangi jumlah reproduksi dasar (\(R_0\)) nilai parameter \(\alpha\) dan \(\beta\) harus dikurangi.

Pada Tabel \ref{tab3} nilai parameter \(\alpha = 1.0563 \times 10^{-8}\) dan \(R_0 = 2.5208\), untuk menjaga stabilitas titik ekuilibrium bebas penyakit (\(E_0\)), maka \(R_0 < 1\). Karena indeks sensitivitas \(\alpha = 1\), agar tidak terjadi penyebaran penyakit tuberkulosis atau endemik, maka diturunkan sebesar 90\% sehingga nilai parameter menjadi \(\alpha = 1.0563 \times 10^{-9}\) dan \(R_0 = 0.2520\). Sementara nilai parameter \(\beta = 0.0010921\) dengan nilai \(R_0 = 2.5208\) dan indeks sensitivitas \(\beta = 0.9327\), agar tidak terjadi penyebaran penyakit tuberkulosis atau endemik, maka diturunkan juga sebesar 90\% sehingga nilai parameter menjadi \(\beta = 0.00010921\) dan \(R_0 = 0.2683\).

Dengan demikian, diperoleh titik kesetimbangan bebas penyakit yang stabil secara asimtotik. Secara medis, untuk mengurangi laju penyebaran penyakit tuberkulosis dari individu yang terinfeksi terhadap individu yang rentan (\(\alpha\)) dan individu yang terinfeksi yang dapat menularkan penyakit (\(\beta\)), dapat dilakukan dengan berbagai tindakan untuk mencegah penularan penyakit tuberkulosis seperti mengisolasi diri, penggunaan masker, kepatuhan minum obat, mengurangi kontak dengan penderita tuberkulosis, dan tindakan pencegahan lainnya. Berikut disajikan gambar mengenai sub populasi exposed, infected, dan removed berdasarkan perubahan parameter \(\alpha\) dan \(\beta\).
\begin{figure}[http]
\centering\includegraphics[width=11cm]{def1}
\caption{Trayektori fase menggunakan nilai parameter pada Tabel \ref{tab2} dengan $\alpha=1.0563 \times 10^{-9}$} \label{gb4}
\end{figure}
\begin{figure}[http]
\centering\includegraphics[width=11cm]{E0b}
\caption{Trayektori fase menggunakan nilai parameter pada Tabel \ref{tab2} dengan $\beta=1.0921 \times 10^{-4}$} \label{gb5}
\end{figure}



\begin{figure}[http]
\centering\includegraphics[width=11cm]{bif}
\caption{Bifurkasi forward pada $R_0 =1$} \label{gb6}
\end{figure}

Gambar \ref{gb4} dan Gambar \ref{gb5} menunjukkan perilaku subpopulasi terpapar (E),  subpopulasi terinfeksi (I), dan subpopulasi pulih (R) menggunakan empat nilai awal yang berbeda dan menurunkan parameter $\alpha$ dan $\beta$ sebesar 90\% akan menuju satu titik ekuilibrium bebas penyakit $E_0=(S,E,I,R)=(8.9174 ,0 ,0 ,0)$.  Hal ini menunjukkan bahwa titik ekuilibrium bebas penyakit $E_0$ bersifat stabil asimtotik. Artinya ketika bilangan reproduksi dasar $R_0<1$ dengan kondisi awal apapun, penyebaran penyakit dari individu terinfeksi sangat rendah sehingga penyakit akan menghilang dari populasi seiring berjalannya waktu. Gambar \ref{gb6} memperlihatkan titik ekuilibrium bebas penyakit $E_0$ dan titik ekuilbrium endemik $E_1$ saling berubah kestabilannya. Perubahan tersebut terjadi saat $R_0=1$ atau bersesuaian dengan $\alpha = \alpha^*=4.1902 \times 10^{-9}$. Titik ekuilibrium bebas penyakit $E_0$ stabil asimtotik lokal jika $R_0 <1$ dan tidak stabil jika $R_0 >1$. Selanjutnya jika $R_0 >1$ terdapat satu titik ekuilibrium endemik $E_1$ yang stabil asimtotik. Pada kondisi terjadi bifurkasi forward ini, pencegahan penyebaran tuberkulosis dapat dilakukan dengan mereduksi $\alpha$ sehingga bernilai kurang dari $\alpha^*$ atau $R_0 <1$.

\section{Kesimpulan}
Pada artikel ini dibahas model SEIR dengan pseudo recovery (sembuh semu) untuk penyakit tuberkulosis. Dari model tersebut diperoleh dua titik ekuilibrium, yaitu titik ekuilibrium bebas penyakit $E_0$ dan titik ekuilibrium endemik $E_1$. Dari analsis kestabilan diperoleh bilangan reproduksi dasar yang diperoleh yaitu $R_0= \frac{\alpha \Lambda \beta \left(\mu +\omega \right)}{\mu \left(\mu + \beta \right)\left(\mu \left(\mu + \delta +\gamma \right)+\omega \left(\mu +\delta \right)\right)}$ sebagai nilai ambang terjadinya endemik. Berdasarkan analisis kestabilan, titik ekuilibrium bebas penyakit stabil asimtotik lokal dan stabil asimtotik global ketika $R_0<1$ dan titik ekuilibrium endemik $E_1$ stabil asimtotik global ketika $R_0>1$.  Berdasarkan analisis sensitivitas, semakin besar nilai parameter $\alpha, \beta,$ dan $\omega$ maka $R_0$ akan semakin besar sehingga meningkatkan kemungkinan terjadi penyebaran penyakit tuberkulosis. Sementara, semakin besar nilai parameter $\mu,\delta,$ dan $\gamma$ maka nilai $R_0$ akan semakin menurun sehingga penyebaran penyakit tuberkulosis juga menurun. Selanjutnya dari analisis bifurkasi, dengan melakukan kontinuasi titik ekuilbrium bebas penyakit terhadap parameter $\alpha$ model ini memiliki bifurkasi forward. Dari hasil tersebut dapat disimpulkan bahwa salah satu upaya pencegahan penyebaran penyakit tuberkulosis adalah dengan menurunkan angka penularan yang dapat dicapai melalui upaya medis seperti menjaga kekebalan tubuh, menggunakan masker, kepatuhan minum obat, dan menghindari kontak dengan penderita tuberkulosis.


\section{Ucapan Terima kasih}
Terima kasih kepada Universitas Ahmad Dahlan yang telah membantu dalam penulisan makalah ini.

\begin{thebibliography}{10}
 \bibitem{1}Irwan, 2019. \emph{Epidemiologi Penyakit Menular}, Edisi ke-3, Absolute Media, Yogyakarta.

\bibitem{2}WHO, 2023.\emph{ Global tuberculosis report 2021}

\bibitem{3}Biro Komunikasi dan Pelayanan Publik, Kementerian Kesehatan RI, \emph{ Deteksi TBC Capai Rekor Tertinggi di Tahun 2022}

\bibitem{4}WHO, 2022. \emph{ Fakta-fakta Utama Tuberkulosis 2022}

\bibitem{6}Aggarwal, A.N., 2019. Quality of life with TB. \emph{Journal of Clinical Tuberculosis and Other Mycobacterial Diseases}, 17. doi: 10.1016/j.jctube.2019.100121.

 \bibitem{7}Ramli, N.A., \& Rahman, H., 2022. Penularan Tuberkulosis Paru dalam Anggota Keluarga di Wilayah Kerja Puskesmas Siko Kota Ternate.

\bibitem{8}Anwari, S.P., 2023. Tren Penyakit TBC Meningkat, Kadinkes Jabar: Masih Banyak Stigma dari Masyarakat.\emph{Portal jabar}.

\bibitem{9}Chakaya, J.M., Harries, A.D., Marks, G.B., 2020. Ending tuberculosis by 2030—Pipe dream or reality?. \emph{International Journal of Infectious Diseases}, 92, S51-S54, doi: 10.1016/j.ijid.2020.02.021.

\bibitem{10}Said, C.S., Side, S., Syam, R., 2020. Model SEIRS Penyebaran Penyakit Tuberkulosis di Kota Makassar.

\bibitem{11}Roudhotillah, D., Chandra, T.D., 2021. Analisis Kestabilan Model Penyebaran Penyakit Tuberkulosis Dengan Menggunakan MSEITR. \emph 15(2), 1858-0629.

\bibitem{12}Qomariyah, N., sutimin, Herdiana, R., Utomo, R.H.S, Permatasari, A.H., 2021. Stability analysis of a tuberculosis epidemic model with nonlinear incidence rate and treatment effects. \emph{ Journal of Physics: Conference Series, IOP Publishing Ltd} , doi: 10.1088/1742-6596/1943/1/012118.

\bibitem{13}Upadhyay, R.K., Pal, A.K., Kumari, S., Roy, P., 2019. Dynamics of an SEIR epidemic model with nonlinear incidence and treatments rates. \emph{ Nonlinear Dyn}, 96(4), 2351-2368, doi: 10.1007/s11071-019-04926-6.

\bibitem{14}Qunia, A.F., Kusnanto, A., Sianturi, P., 2021. The Impact of Relapse Rate on Deterministic Epidemiological Models with Pseudo-recovery. \emph{ Journal of Physics: Conference Series, IOP Publishing Ltd}, doi: 10.1088/1742-6596/1863/1/012004.

\bibitem{15}Perko, 2001, \emph{Texts in Applied Mathematics 7: Differential Equations and Dynamical System}, Third Edition, Springer, New York.

\bibitem{16}Murray, J.D., 2002, \emph{Mathematical Biology: I. An Introduction}, Third Edition, Springer 2002.

\bibitem{17}Yeketi, A. A., \& Othman, W. A. M., Musa, R., 2021. A Mathematical Model of the Tuberculosis Epidemic.\emph{Acta Biotheor}, 69(3), 225-255. doi: 10.1007/s10441-020-09406-8.

\bibitem{18}Sulistiyowati, A., Abadi, 2023. Analisis Kestabilan Model Penyebaran Tuberkulosis Dengan MDR-TB dan Pengaruh Vaksinasi. \emph{Jurnal Ilmiah Matematika}, 11(2).

\bibitem{19}Ihsan, H., Side, M., Pagga, J., 2021. Pemodelan Matematika SEIRS Pada Penyebaran Penyakit Malaria di Kabupaten Mimika.

\bibitem{20}Panigoro, H.S., Savitri, D., 2020. Bifurkasi Hopf pada Model Lotka-Volterra Orde-Fraksional dengan Efek Allee Aditif pada Predator. \emph{Jambura Journal of Biomathematics (JJBM)}, 1(1), 16-24, doi: 10.34312/jjbm.v1i1.6908.

\bibitem{21}Maligranda, L., 2012. The AM-GM Inequality is Equivalent to the Bernoulli Inequality.\emph{ Math Intelligencer}, 34, 1-2, https://doi.org/10.1007/s00283-011-9266-8.

\bibitem{22}BPS, 2023. \emph{Jumlah Penduduk Menurut Kabupaten/Kota (Jiwa)}, Badan Pusat Statistika (BPS).

\bibitem{23}Lestari, R.,\emph{Percepatan Eliminasi TBC Difokuskan Pada Peningkatan Penemuan Kasus}. https://diskes.jabarprov.go.id/informasipublik/detail-berita. Diakses pada 25 Maret 2024.

 \bibitem{24}Castillo-Chavez, C., Song, B., 2020. Dynamical models of tuberculosis and their applications. \emph{Mathematical Biosciences and Engineering} 1(2): 361–404. doi:DOI:10.3934/mbe.2004.1.361.
\bibitem{25}Chitnis, N., Hayman, J. M.,  Cushing, J. M., 2008. Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. \emph{ Bull Math Biol}, 70(5), 1272-1296, doi: 10.1007/s11538-008-9299-0.



\end{thebibliography}
\end{document}
