\documentclass{template-jurnal} %Bagian ini diedit oleh editor Jurnal Matematika UNAND

\hyphenation{di-tulis-kan di-terang-kan}
%Perhatikan aturan penulisan dan ukuran huruf yang digunakan
\begin{document}

\markboth{Nisardi, dkk} %Jika lebih dari dua penulis, tuliskan sebagai Nama Penulis Pertama dkk.
{A Deterministic Mathematical Model of Meningitis Transmission Dynamics with Vaccination and Screening}

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

\title{A Deterministic Mathematical Model of Meningitis Transmission Dynamics with Vaccination and Screening}

%\author{Muhammad Rifki Nisardi$^{a}$\footnote{penulis korespondensi}, Muh Ikhsan Amar$^{a}$, Muhammad Fadhil Nurahmad\\ }

%\address{ Institut Teknologi Bacharuddin Jusuf Habibie, Parepare, Indonesia,\\
\\
\\
%email : \email{muhammadrifkinisardi@ith.ac.id,\\ ikhsan.amar93@ith.ac.id, \\
%nmuhfadhil@gmai.com}}

% Jika instansi semua penulis sama, cukup tuliskan satu kali saja, tanpa menggunakan $^{a}$

\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
%Abstrak tidak melebihi 250 kata. Tuliskan intisari dari makalah anda.

%\keywords{Tuliskan maksimal tiga buah kata kunci, tuliskan sesuai urutan abjad}
%\bigskip

\textbf{Abstract}. % Dalam bahasa Inggris
\textit{This study aims to examine the mathematical model of meningitis transmission as a deterministic model. The model includes five compartments: susceptible (S), carrier (C), infected (I), treatment (T), and recovered (R). We also consider vaccination and screening as interventions in disease transmission. In this work, we obtained two equilibrium points: disease-free equilibrium point and endemic equilibrium point. The next generation matrix is employed to compute the basic reproduction numbers ($R_0$). We also analyzed the sensitivity of parameters concerning $R_0$. If $R_0 < 1$, then the disease-free equilibrium point exists and is locally stable, whereas the endemic equilibrium point exists when $R_0 > 1$ and is locally stable if the Routh-Hurwitz criterion is satisfied. We use the Runge-Kutta 4th order method to confirm the stability properties of the equilibrium points and also show that vaccination and screening affect the transmission dynamics of Meningitis.}

\noindent \textit{Keywords}: Basic Reproduction Number, Meningitis Modelling Disease, Routh-Hurwitz Criterion
\end{abstract}


\section{Introduction}
Mathematical models are utilized as tools to depict the behavior of real-world problems under various assumptions, particularly in epidemiological issues. Numerous studies have modeled disease transmission through mathematical frameworks \cite{ref1, ref2, ref3, ref4, ref5, ref6, ref7, ref8}. These studies aim to provide a theoretical understanding of the dynamics of specific disease spread, thereby serving as a reference to assist in decision-making for addressing epidemiological problems within a population. One disease with a high risk of transmission is Meningitis.


Meningitis is an infectious disease caused by the bacterium \textit{Neisseria meningitidis} (meningococcus) \cite{ref9}. Bacterial meningitis infection can occur through the transfer of bacteria contained in respiratory secretions from an infected individual to a susceptible individual via airborne transmission or the use of contaminated personal items. Upon entering the body, \textit{Neisseria meningitidis} infiltrates the bloodstream and infects the meninges \cite{ref10, ref11, ref12, ref13}. The meninges are thin layers or membranes that protect the brain and spinal cord \cite{ref14}. Infected individuals may be asymptomatic or symptomatic, showing symptoms such as high fever, headache, stiff neck, vomiting, and skin rash. If not promptly treated, meningitis can cause swelling of the fluid around the brain and spinal cord, potentially leading to disability or death \cite{ref13, ref15}.


The spread of Meningococcal Meningitis is global, with 6,469 cases reported in 2023, including 570 confirmed cases and 420 deaths due to the infection. According to data from the Indonesian Ministry of Health, one in five Meningococcal Meningitis cases results in long-term disabilities, such as limb loss, deafness, neurological issues, and brain damage. The Case Fatality Rate (CFR) for Meningococcal Meningitis ranges from 5-15\%, depending on clinical symptoms \cite{ref16}. In endemic regions, there are approximately 10 cases per 100,000 population annually. In 2016, Indonesia had the highest incidence and mortality rate of Meningitis in Southeast Asia \cite{ref9}.

This study focuses on a deterministic mathematical model of the dynamics of Meningitis transmission. The model divides the population into six compartments based on different clinical conditions. Furthermore, this study’s model considers control factors such as vaccination and screening, which can influence disease transmission dynamics. The results of this study are presented in an article comprising several sections. 

\section{Literature Review}
Several previous studies are pertinent. The first study \cite{ref17} investigates a mathematical model of the dynamics of Meningitis transmission, involving two control measures: vaccination and treatment. The researchers developed a mathematical model with four compartments: susceptible, carrier, infected, and recovered. In addition to model analysis, this study addresses the optimal control problem with vaccination and treatment rates as control parameters.

Next, \cite{ref10} examines a Meningitis model focusing on cases in Northern Nigeria. The researchers constructed an ordinary differential equation model consisting of six compartments: susceptible, carrier of infection 1, carrier of infection 2, infected with strain 1, infected with strain 2, and recovered. The model identifies three equilibrium points: the disease-free equilibrium ($E_0$), the endemic equilibrium for strain 1 ($E_1$), and the endemic equilibrium for strain 2 ($E_2$).

Another study \cite{ref18} focuses on a Meningitis transmission model involving influenza co-infection and considers the quarantine of individuals infected with influenza. The researchers divided the population into six sub-populations: susceptible, influenza infectives only, meningitis infectives only, influenza-meningitis co-infectives, quarantine of influenza infectives, and recovered. The study's findings indicate that increasing the quarantine rate of individuals infected with influenza and the recovery rate will reduce the number of individuals infected with influenza, meningitis, or both.

Additionally, \cite{ref19} examines a mathematical model of the spread of co-infections of meningitis and pneumonia, incorporating interventions that can help reduce the number of infections within a population. The researchers developed a deterministic ordinary differential equation model and divided the population into seven sub-populations: susceptible, individuals with pneumonia infection, individuals with meningitis infection, individuals with co-infections of meningitis and pneumonia, co-infected individuals receiving treatment, vaccinated individuals, and recovered individuals. The study's results show that vaccination, reducing contact with infected individuals, and treatment have a significant impact on reducing the mortality rate from these diseases.


\section{Model Formulation}
In this study, the Meningitis transmission model is divided into five subpopulations based on different clinical conditions: susceptible $S(t)$ representing susceptible individuals; carrier $C(t)$ representing infected individuals without clinical symptoms but capable of spreading the disease; infected $I(t)$ representing individuals with clinical symptoms; treatment $T(t)$ representing individuals receiving treatment; and recovery $R(t)$ representing individuals who have recovered from meningitis.

\begin{eqnarray}
    \frac{dS}{dt} &=& \Lambda - \mu S - \frac{\beta(1 - \epsilon v)(\eta C + \eta_1 I) S}{N} - v \epsilon S + \theta R , \nonumber\\
    \frac{dC}{dt} &=& \frac{\beta(1 - \epsilon v)(\eta C + \eta_1 I)}{N} S - (\mu + \omega + \gamma + \sigma) C , \nonumber\\
    \frac{dI}{dt} &=& \sigma C - (\delta + \mu + \gamma + \phi) I , \nonumber\\ 
    \frac{dT}{dt} &=& \omega C + \phi I - (\delta + \mu + \gamma_1) T , \nonumber \label{eq:basicmodel}\\ 
    \frac{dR}{dt} &=& \gamma_1 T + \gamma I + \gamma C - (\mu + \theta) R + v \epsilon S, \nonumber \\ 
\end{eqnarray}



In Model (\ref{eq:basicmodel}), the individuals in the $S$ compartment increase due to natural births at a rate $\Lambda$, and individuals who lose immunity at a rate $\theta$. They decrease due to vaccination at a rate $v$ with efficacy $\epsilon$, natural death at a rate $\mu$, and interaction between susceptible individuals and infected or carrier individuals with a contact transmission rate $\beta$, causing susceptible individuals to move to the $C$ compartment.

The $C$ compartment increases due to susceptible individuals infected by carriers or infected individuals with proportions $\eta$ and $\eta_1$ respectively. This compartment decreases due to natural death, transition to the infected compartment, natural recovery, and early tracing for treatment.

Next, the $I$ compartment increases due to the transition from the $C$ compartment at a rate $\sigma$. This compartment decreases due to natural death, meningitis-induced death, natural recovery, and transition to the treatment compartment. The treatment compartment increases due to individuals traced from the $C$ compartment at a rate $\omega$ and infected individuals receiving treatment at a rate $\phi$. This compartment decreases due to natural death, meningitis-induced death, and recovery at a rate $\gamma_1$. Finally, the recovery compartment increases due to individuals recovering from meningitis and those receiving vaccination. This compartment decreases due to natural death and loss of immunity at a rate $\theta$. The interactions among compartments are illustrated in Figure 1.
\begin{figure}[htbp]
    \centering
    \includegraphics[width=0.8\linewidth]{Model Meningitis.png}
    \caption{Flowchart of Model}
    \label{fig:flow-model}
\end{figure}


Each parameter in Model (1) is assumed to have a positive value.

\begin{table}[h]
\centering
\begin{tabular}{|c|l|}
\hline
\textbf{Parameter} & \textbf{Description} \\ \hline
$\Lambda$ & Recruitment rate \\ \hline
$\mu$ & Natural death rate \\ \hline
$\beta$ & Effective contact rate \\ \hline
$\eta$ & Per capita infection rate by carriers \\ \hline
$\eta_1$ & Per capita infection rate by ill individuals \\ \hline
$v$ & Vaccination rate against meningitis \\ \hline
$\theta$ & Loss of immunity \\ \hline
$\omega$ & Tracing rate \\ \hline
$\sigma$ & Rate of progression from C to I \\ \hline
$\varphi$ & Treatment rate \\ \hline
$\delta$ & Disease-induced mortality \\ \hline
$\gamma$ & Natural recovery rate from disease \\ \hline
$\gamma_1$ & Recovery rate with treatment \\ \hline
\end{tabular}
\caption{Parameter Descriptions}
\end{table}

\subsection{Analysis of Positive (Non-Negativity) and Boundedness of Solutions}
The system (1) involves interactions between populations, thus the solution of the system must be positive and bounded. For system (1), a feasible region $\omega \subseteq \mathbb{R}_+^5$ is defined such that 
\[
\omega = \{(S,C,I,T,R) \in \mathbb{R}_+^5 : N(t) \leq \frac{\Lambda}{\mu}\}
\]

\textbf{Theorem 1.} The solution of system (1) within $\Omega \subseteq \mathbb{R}_+^5$ is non-negative and bounded.

\textbf{Proof:} From the first equation of system (1),
\[
\frac{dS}{dt} = \lambda - \mu S - \lambda S - \epsilon v S + \theta R \geq -(\mu + \lambda + \epsilon v)S \tag{3.2}
\]
Thus,
\[
\frac{1}{S} \frac{dS}{dt} \geq -(\mu + \lambda + \epsilon v) \tag{3.3}
\]
Integrating equation (3.3) yields
\[
\ln S \geq -(\mu + \lambda + v)t + K, \tag{3.4}
\]
where \( K \) is the integration constant. For the initial condition \( S(0) = S_0 \), we obtain
\[
S \geq S_0 e^{-(\mu + \lambda + v)t} \geq 0. \tag{3.5}
\]
Similarly, it can be shown that the other equations in system \ref{eq:basicmodel} are also positive for all \( t > 0 \). Next, we will show that the solution of system (1) is bounded. It is known that \( N(t) = S(t) + C(t) + I(t) + T(t) + R(t) \). Based on system (1), we get
\[
\frac{dN}{dt} = \frac{dS}{dt} + \frac{dC}{dt} + \frac{dI}{dt} + \frac{dT}{dt} + \frac{dR}{dt}
\]
\[
\frac{dN}{dt} = \Lambda - \mu(S + C + I + T + R) - \delta(I + T)
\]
\[
\frac{dN}{dt} = \Lambda - \mu N - \delta(I + T) \leq \Lambda - \mu N \tag{3.6}
\]
Since \( N(t) \geq 0 \), then \(\frac{dN}{dt} \geq 0\), so we obtain \( 0 \leq \Lambda - \mu N \) or \( N(t) \leq \frac{\Lambda}{\mu} \). Therefore, it is proven that the solution of system \ref{eq:basicmodel} is bounded. 

\subsection{Normalization of the Model}
Model \ref{eq:basicmodel} involves interactions among individuals in a population. To simplify the calculations and simulations, each number of individuals in the compartments is expressed as a proportion of the total population as follows:
\[
x_1 = \frac{S}{N}, \quad x_2 = \frac{C}{N}, \quad x_3 = \frac{I}{N}, \quad x_4 = \frac{T}{N}, \quad x_5 = \frac{R}{N} \tag{3.7}
\]
By substituting Equation (3.7) into Model \ref{eq:basicmodel}, we obtain:
\begin{align*}
\frac{dx_1}{dt} &= \mu - \mu x_1 - \beta (1 - \epsilon v) (\eta x_2 + \eta_1 x_3) x_1 - v \epsilon x_1 + \theta x_5 \\
\frac{dx_2}{dt} &= \beta (1 - \epsilon v) (\eta x_2 + \eta_1 x_3) x_1 - (\mu + \omega + \gamma + \sigma) x_2 \\
\frac{dx_3}{dt} &= \sigma x_2 - (\delta + \mu + \gamma + \varphi) x_3 \\
\frac{dx_4}{dt} &= \omega x_2 + \varphi x_3 - (\delta + \mu + \gamma_1) x_4 \\
\frac{dx_5}{dt} &= \gamma_1 x_4 + \gamma x_3 + \gamma x_2 - (\mu + \theta) x_5 + v \epsilon x_1 \tag{3.8}
\end{align*}

\section{Mathematical Analysis and Numerical Result}
\subsection{Equilibrium Points}

In the case of the meningitis transmission model, there are two equilibrium points in the system: the disease-free equilibrium point and the endemic equilibrium point. The disease-free equilibrium point occurs when \( C = I = 0 \), which means that there are no infected individuals in the system spreading the disease to other individuals in the population. By assuming \( C = I = 0 \) and substituting these values into System (3.8), the disease-free equilibrium point is obtained as follows:
\[
E^0 = \left( x_1^0, x_2^0, x_3^0, x_4^0, x_5^0 \right) = \left( \frac{\mu + \theta}{\mu + \theta + \epsilon v}, 0, 0, 0, \frac{\epsilon v}{\mu + \theta + \epsilon v} \right)
\]

Next, the condition where there are infected individuals who can transmit meningitis to others in the population is referred to as the endemic condition. In other words, the endemic condition occurs when \( I \neq 0 \). If the system's solution tends to move towards the endemic equilibrium point, it means that meningitis will continue to exist in the population. From the analysis, the endemic equilibrium point for model (3.8) is obtained as \( E^1 = \left( x_1^*, x_2^*, x_3^*, x_4^*, x_5^* \right) \). The values of \( S^*, C^*, T^*, R^* \) are expressed in Equations (9)–(13) and depend on \( I^* \) as follows:

\[
x_1^* = \frac{\mu \sigma m_2 (\mu + \theta) + \theta (\gamma_1 (\omega m_1 + \sigma \varphi) + \gamma m_2 (\sigma + m_1)) x_3^*}{\mu \sigma (\mu + \theta + \epsilon v) + \beta (1 - \epsilon v)(\mu + \theta)(\eta m_1 + \eta_1 \sigma) x_3^*} \tag{3.9}
\]

\[
x_2^* = \frac{m_1}{\sigma} x_3^* \tag{3.10}
\]

\[
x_3^* = \frac{\mu \sigma \left\{ m_1 m_3 (\mu + \theta + \epsilon v) - \beta m_2 (1 - \epsilon v)(\mu + \theta)[\sigma \eta_1 + \eta m_1 ] \right\}}{\beta (1 - \epsilon v)[\sigma \eta_1 + \eta m_1 ]\left\{ \theta (\gamma_1 (\omega m_1 + \sigma \varphi) + \gamma \sigma m_2 + \gamma m_1 m_2 ) - (\mu + \theta) m_1 m_3 \right\}} \tag{3.11}
\]

\[
x_4^* = \frac{\omega m_1 + \varphi \sigma}{\sigma m_2} x_3^* \tag{3.12}
\]

\[
x_5^* = \frac{P_1 {x_3^*}^2 + P_2 x_3^* + P_3}{\sigma m_2 (\mu + \theta)[\mu \sigma (\mu + \theta + \epsilon v) + \beta (1 - \epsilon v)(\mu + \theta)(\eta m_1 + \eta_1 \sigma) x_3^*]} \tag{3.13}
\]

where\\


\(m_1 = \delta + \mu + \gamma + \varphi \)\\
\(m_2 = \delta + \mu + \gamma_1\)\\
\(m_3 = \mu + \omega + \gamma + \sigma\)\\
\(P_1 = (\gamma_1 (\omega m_1 + \sigma \varphi) + \gamma m_2 (\sigma + m_1))(\beta (1 - \epsilon v)(\mu + \theta)(\eta m_1 + \eta_1 \sigma))\)\\
\(P_2 = \sigma (\gamma_1 (\omega m_1 + \sigma \varphi) + \gamma m_2 (\sigma + m_1))(\mu (\mu + \theta + \epsilon v) + \theta \epsilon v m_2 )\)\\
\(P_3 = (m_2 )^2 \sigma^2 \mu \epsilon v (\mu + \theta)\)\\
\subsection{Basic Reproduction Number}

The basic reproduction number will be determined using the Next Generation Matrix (NGM) method \cite{ref20}. Let $F_i(x)$ represent the rate of new infections in compartment $i$ and $V_i(x)$ represent the rate of transfer of individuals in compartment $i$. Thus, $F_i(x)$ and $V_i(x)$ are given by:

\begin{equation}
F_i(x) = \begin{bmatrix}
\beta (1 - \epsilon v) (\eta x_2 + \eta_1 x_3) x_1 \\
0 \\
0 
\end{bmatrix},
\end{equation}

\begin{equation}
V_i(x) = \begin{bmatrix}
(\mu + \omega + \gamma + \sigma) x_2 \\
(\delta + \mu + \gamma + \phi) x_3 - \sigma x_2 \\
(\delta + \mu + \gamma_1) x_4 - \omega x_2 - \phi x_3 
\end{bmatrix}.
\end{equation}

From Equations (4.1) and (4.2), we obtain the matrices $F$ and $V$, which are then evaluated at the disease-free equilibrium point $E_0$ as follows:

\begin{equation}
F = \frac{\partial F_i(x)}{\partial (x_2, x_3, x_4)} = \begin{bmatrix}
\beta (1 - \epsilon v) \eta x_1^0 & \beta (1 - \epsilon v) \eta_1 x_1^0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix},
\end{equation}

\begin{equation}
V = \frac{\partial V_i(x)}{\partial (x_2, x_3, x_4)} = \begin{bmatrix}
\mu + \omega + \gamma + \sigma & 0 & 0 \\
-\sigma & \delta + \mu + \gamma + \phi & 0 \\
-\omega & -\phi & \delta + \mu + \gamma_1
\end{bmatrix} = \begin{bmatrix}
k_1 & 0 & 0 \\
-\sigma & k_2 & 0 \\
-\omega & -\phi & k_3
\end{bmatrix}.
\end{equation}

The inverse of matrix $V$ is:

\begin{equation}
V^{-1} = \begin{bmatrix}
1/k_1 & 0 & 0 \\
\sigma / (k_1 k_2) & 1/k_2 & 0 \\
(\omega k_2 + \phi \sigma) / (k_1 k_2 k_3) & \phi / (k_2 k_3) & 1/k_3
\end{bmatrix}.
\end{equation}

Next, using the matrices $F$ and $V^{-1}$, we obtain the matrix $FV^{-1}$:

\begin{equation}
FV^{-1} = \begin{bmatrix}
\frac{\beta (1 - \epsilon v) \eta x_1^0}{k_1} + \frac{\beta (1 - \epsilon v) \eta_1 x_1^0 \sigma}{k_1 k_2} & \frac{\beta (1 - \epsilon v) \eta_1 x_1^0}{k_2} & 0 \\
0 & 0 & 0 \\
0 & 0 & 0
\end{bmatrix}.
\end{equation}

Since the basic reproduction number $R_0$ is the spectral radius of $FV^{-1}$, or in other words, $R_0 = \max(\lambda_i)$, we obtain the basic reproduction number:

\begin{equation}
R_0 = \frac{\beta (1 - \epsilon v) (\mu + \theta) (\eta (\delta + \mu + \gamma + \phi) + \eta_1 \sigma)}{(\mu + \omega + \gamma + \sigma)(\delta + \mu + \gamma + \phi)(\mu + \epsilon v + \theta)}.
\end{equation}

Equation (4.7) is then plotted to illustrate the parameter combinations' effect on $R_0$. This simulation is conducted by varying two parameters and fixing the others. The parameters include combinations of $\beta$, $\epsilon$, $v$, and $\omega$. The simulation results are provided in Figure 2.
\begin{figure}[htbp]
    \centering
    \includegraphics[width=1\linewidth]{Profil R0.png}
    \caption{Basic Reproduction Number Profile}
    \label{fig:R0 profile}
\end{figure}

\subsection{Sensitivity Analysis}

The sensitivity analysis of the basic reproduction number is conducted to determine the influence of each parameter on the basic reproduction number $R_0$. This analysis focuses on the parameters present in the basic reproduction number $R_0$. The objective is to observe the effect of parameter variations on the basic reproduction number $R_0$. The normalized sensitivity index is obtained by normalizing the differential of the variable $R_0$ with respect to parameter $p$, defined as follows \cite{ref21}:
\begin{align}
    C_p^{R_0} = \left(\frac{\partial R_0}{\partial p}\right) \times \frac{p}{R_0}
\end{align}



In this case, the variable to be measured is the basic reproduction number of the Meningitis transmission model with respect to each parameter. The results of the parameter sensitivity analysis are as follows:

\begin{eqnarray}
C_\beta^{R_0} &=& 1, \nonumber\\
C_v^{R_0} &=& -\frac{\mu\epsilon + \theta\epsilon + \epsilon v}{(1 - \epsilon v)(\epsilon v + \mu + \theta)}, \nonumber\\
C_{\gamma_1}^{R_0} &=& 0, \nonumber\\
C_\theta^{R_0} &=& \frac{\theta \epsilon v}{(\epsilon v + \mu + \theta)(\mu + \theta)}, \nonumber\\
C_\eta^{R_0} &=& \frac{\eta (\delta + \mu + \gamma + \phi)}{\eta (\delta + \mu + \gamma + \phi) + \eta_1 \sigma}, \nonumber\\
C_\delta^{R_0} &=& -\frac{\eta_1 \sigma \delta}{(\delta + \mu + \gamma + \phi)(\eta (\delta + \mu + \gamma + \phi) + \eta_1 \sigma)}, \nonumber\\
C_\gamma^{R_0} &=& -\frac{(\delta + \mu + \gamma + \phi)^2 \eta + 2 (\gamma + \mu + \frac{1}{2}\phi + \frac{1}{2}\delta + \frac{1}{2}\omega + \frac{1}{2}\sigma) \sigma \eta_1 \gamma}{(\gamma + \mu + \omega + \sigma)(\delta + \mu + \gamma + \phi)(\eta (\delta + \mu + \gamma + \phi) + \eta_1 \sigma)}, \nonumber\\
C_\phi^{R_0} &=& -\frac{\phi \eta_1 \sigma}{(\delta + \mu + \gamma + \phi)(\eta (\delta + \mu + \gamma + \phi) + \eta_1 \sigma)}, \nonumber\\
C_{\eta_1}^{R_0} &=& \frac{\sigma \eta_1}{\eta (\delta + \mu + \gamma + \phi) + \eta_1 \sigma}, \nonumber\\
C_\sigma^{R_0} &=& -\frac{(\eta (\delta + \mu + \gamma + \phi) - \eta_1 (\gamma + \mu + \omega)) \sigma}{(\gamma + \mu + \omega + \sigma)(\eta (\delta + \mu + \gamma + \phi) + \eta_1 \sigma)}, \nonumber\\
C_\omega^{R_0} &=& -\frac{\omega}{(\gamma + \mu + \omega + \sigma)}, \nonumber\\
C_\epsilon^{R_0} &=& -\frac{\mu\epsilon + \theta\epsilon + \epsilon v}{(\mu + \theta + \epsilon v)(1 - \epsilon v)}, \nonumber\\
\end{eqnarray}

\begin{table}[h!]
\centering
\begin{tabular}{|c|c|c|}
\hline
Parameter & Description & Sensitivity Index \\
\hline
$\beta$ & Effective Contact Rate & 1 \\
$v$ & Vaccination Rate against Meningitis & -2.0763 \\
$\theta$ & Loss of Immunity & 0.6796 \\
$\omega$ & Tracing Rate & -0.8561 \\
$\sigma$ & Rate of Progression from C to I & 0.0831 \\
$\phi$ & Treatment Rate & -0.0858 \\
$\delta$ & Disease-Induced Mortality & -0.0028 \\
$\gamma$ & Natural Recovery Rate from Disease & -0.1257 \\
$\gamma_1$ & Recovery Rate with Treatment & 0 \\
$\epsilon$ & Effectiveness of Vaccination & -2.0763 \\
$\eta$ & Per Capita Infection Rate by Carriers & 0.8483 \\
$\eta_1$ & Per Capita Infection Rate by Ill Individuals & 0.1516 \\
\hline
\end{tabular}
\caption{Sensitivity Index $R_0$ to the parameters}
\end{table}
\begin{figure}
    \centering
    \includegraphics[width=0.75\linewidth]{Sens_Analysis.png}
    \caption{Sensitivity Analysis Index}
    \label{fig:SensIndex}
\end{figure}



\subsection{Local Stability Analysis of Equilibrium Points}
If the system (3.8) is linearized, a Jacobian matrix is obtained as follows:
\begin{equation}
J = \begin{pmatrix}
-(\mu + \lambda + \varepsilon v) & \beta (1 - \varepsilon v) \eta x_1 & \beta (1 - \varepsilon v) \eta_1 x_1 & 0 & \theta \\
\lambda & \beta (1 - \varepsilon v) \eta x_1 - (\mu + \omega + \gamma + \sigma) & \beta (1 - \varepsilon v) \eta_1 x_1 & 0 & 0 \\
0 & \sigma & -(\delta + \mu + \gamma + \phi) & 0 & 0 \\
0 & \omega & \phi & -(\delta + \mu + \gamma_1) & 0 \\
-\varepsilon v & \gamma & \gamma & \gamma & -(\mu + \theta) 
\end{pmatrix}
\end{equation}

\begin{itemize}
    \item {Local Stability of the Disease-Free Equilibrium Point}
\end{itemize}

Lemma: The disease-free equilibrium point \(E_0\) in system (3.8) is locally asymptotically stable if \(R_0 < 1\) and unstable if \(R_0 > 1\).

Proof:
\begin{equation}
J(E_0) = \begin{pmatrix}
-(\mu + \varepsilon v) & \beta (1 - \varepsilon v) \eta x_1^0 & \beta (1 - \varepsilon v) \eta_1 x_1^0 & 0 & \theta \\
0 & \beta (1 - \varepsilon v) \eta x_1^0 - (\mu + \omega + \gamma + \sigma) & \beta (1 - \varepsilon v) \eta_1 x_1^0 & 0 & 0 \\
0 & \sigma & -(\delta + \mu + \gamma + \phi) & 0 & 0 \\
0 & \omega & \phi & -(\delta + \mu + \gamma_1) & 0 \\
\varepsilon v & \gamma & \gamma & \gamma_1 & -(\mu + \theta) 
\end{pmatrix}
\end{equation}

To simplify notation, Equation (4.11) is written as:
\begin{equation}
J(E_0) = \begin{pmatrix}
-k_1 & h_1 & h_2 & 0 & \theta \\
0 & h_1 - k_2 & h_2 & 0 & 0 \\
0 & \sigma & -k_3 & 0 & 0 \\
0 & \omega & \phi & -k_4 & 0 \\
\varepsilon v & \gamma & \gamma & \gamma_1 & -k_5 
\end{pmatrix}
\end{equation}
where:
\[
k_1 = \mu + \varepsilon v; \quad k_2 = \mu + \omega + \gamma + \sigma; \quad k_3 = \delta + \mu + \gamma + \phi; \quad k_4 = \delta + \mu + \gamma_1; \quad k_5 = \mu + \theta
\]
\[
h_1 = \beta (1 - \varepsilon v) \eta x_1^0; \quad h_2 = \beta (1 - \varepsilon v) \eta_1 x_1^0
\]

From the Jacobian matrix in Equation (4.12), the characteristic equation is obtained:
\begin{equation}
\det(\lambda I - J(E_0)) = 0
\end{equation}
\begin{equation}
(\lambda - k_4)(\lambda^2 + (k_1 + k_5) \lambda + k_1 k_5 + \theta \varepsilon v)(\lambda^2 + (k_2 + k_3 - h_1) \lambda + k_2 k_3 - h_1 k_3 - h_2 \sigma) = 0
\end{equation}

From the characteristic equation (4.14), the eigenvalues are obtained as follows:
\begin{equation}
\lambda_1 = -k_4
\end{equation}
Since every parameter used in Equation (4.15) is positive, this results in \(\lambda_1\) being negative. The stability of the disease-free equilibrium point depends on the eigenvalues \(\lambda_2, \lambda_3, \lambda_4,\) and \(\lambda_5\) obtained from the equations:
\begin{equation}
A(\lambda) = \lambda^2 + A_1 \lambda + A_2 = 0
\end{equation}
\begin{equation}
B(\lambda) = \lambda^2 + B_1 \lambda + B_2 = 0
\end{equation}

Based on the Routh-Hurwitz criteria, Equation (4.16) has eigenvalues with negative real parts if and only if it satisfies:
\[ D_1 > 0, \; D_2 > 0 \]

For the first condition:
\[ D_1 = A_1 \]
\[ D_1 = 2\mu + \varepsilon v + \theta \]
Since every parameter is positive, it is clear that \(D_1 > 0\).

For the second condition:
\[ D_2 = \begin{vmatrix} A_1 & 0 \\ 1 & A_2 \end{vmatrix} \]

\[ D_2 = A_1 A_2 > 0 \]
Since \(A_1, A_2 > 0\), \(D_2 > 0\).

Similarly, for Equation (4.17), the eigenvalues have negative real parts if and only if it satisfies:
\[ D_1 > 0, \; D_2 > 0 \]

For the first condition:
\[ D_1 = B_1 \]
\[ D_1 = 2\mu + \omega + 2\gamma + \sigma + \delta + \phi - \beta (1 - \varepsilon v) \eta x_1^0 \]
Since every parameter is positive, \(D_1 > 0\) if and only if:
\[ 2\mu + \omega + 2\gamma + \sigma + \delta + \phi > \beta (1 - \varepsilon v) \eta x_1^0 \]

For the second condition:
\[ D_2 = \begin{vmatrix} B_1 & 0 \\ 1 & B_2 \end{vmatrix} \]

\[ D_2 = B_1 B_2 > 0 \]
If the condition \(B_1 > 0\) is met, then \(D_2 > 0\) can be satisfied if and only if \(B_2 > 0\) or
\[
(\mu + \omega + \gamma + \sigma)(\delta + \mu + \gamma + \phi) > \beta (1 - \varepsilon v) \eta x_1^0 (\delta + \mu + \gamma + \phi) + \beta (1 - \varepsilon v) \eta_1 x_1^0 \sigma
\]
\[
(\mu + \omega + \gamma + \sigma)(\delta + \mu + \gamma + \phi)(\mu + \varepsilon v + \theta) > \beta (1 - \varepsilon v)(\mu + \theta)(\eta (\delta + \mu + \gamma + \phi) + \eta_1 \sigma)
\]
\[
\frac{\beta (1 - \varepsilon v)(\mu + \theta)(\eta (\delta + \mu + \gamma + \phi) + \eta_1 \sigma)}{(\mu + \omega + \gamma + \sigma)(\delta + \mu + \gamma + \phi)(\mu + \varepsilon v + \theta)} < \frac{(\mu + \omega + \gamma + \sigma)(\delta + \mu + \gamma + \phi)(\mu + \varepsilon v + \theta)}{(\mu + \omega + \gamma + \sigma)(\delta + \mu + \gamma + \phi)(\mu + \varepsilon v + \theta)} = R_0 < 1
\]
Thus, the condition in Equation (4.17) can be met if \(R_0 < 1\).

\begin{itemize}
    \item {Stability of the Endemic Equilibrium Point}
\end{itemize}

\begin{equation}
J(E^1) = \begin{pmatrix}
-(\mu + \lambda + \varepsilon v) & \beta (1 - \varepsilon v) \eta x_1^* & \beta (1 - \varepsilon v) \eta_1 x_1^* & 0 & \theta \\
\lambda & \beta (1 - \varepsilon v) \eta x_1^* - (\mu + \omega + \gamma + \sigma) & \beta (1 - \varepsilon v) \eta_1 x_1^* & 0 & 0 \\
0 & \sigma & -(\delta + \mu + \gamma + \phi) & 0 & 0 \\
0 & \omega & \phi & -(\delta + \mu + \gamma_1) & 0 \\
\varepsilon v & \gamma & \gamma & \gamma & -(\mu + \theta) 
\end{pmatrix}
\end{equation}

Let:
\[
m_1 = \mu + \varepsilon v + \lambda; \quad m_2 = \mu + \omega + \gamma + \sigma; \quad m_3 = \delta + \mu + \gamma + \phi; \quad m_4 = \delta + \mu + \gamma_1; \quad m_5 = \mu + \theta
\]
\[
g_1 = \beta (1 - \varepsilon v) \eta x_1^*; \quad g_2 = \beta (1 - \varepsilon v) \eta_1 x_1^*
\]

To simplify notation, substituting these into matrix (4.18), we obtain:
\begin{equation}
J(E^1) = \begin{pmatrix}
-m_1 & g_1 & g_2 & 0 & \theta \\
\lambda & g_1 - m_2 & g_2 & 0 & 0 \\
0 & \sigma & -m_3 & 0 & 0 \\
0 & \omega & \phi & -m_4 & 0 \\
-\varepsilon v & \gamma & \gamma & \gamma & -m_5 
\end{pmatrix}
\end{equation}

The characteristic equation is obtained as follows:
\begin{equation}
P(\xi) = \xi^5 + a_1 \xi^4 + a_2 \xi^3 + a_3 \xi^2 + a_4 \xi + a_5
\end{equation}
where:\\
\begin{itemize}
    \item \(
     a_1 = m_1 + m_2 + m_3 + m_4 + m_5 - g_1
     \)
    \item \(a_2 = m_1 m_2 + m_1 m_3 + m_1 m_4 + m_1 m_5 + m_2 m_3 + m_2 m_4 + m_2 m_5 + m_3 m_4 + m_3 m_5 + m_4 m_5 + \theta \varepsilon v - g_1 m_1 - g_1 m_3 - g_1 m_4 - g_1 m_1 - g_1 \lambda - g_2 \sigma \)
    
    \item \(a_3 = m_1 m_2 m_3 + m_1 m_2 m_4 + m_1 m_2 m_5 + m_1 m_3 m_4 + m_1 m_3 m_5 + m_1 m_4 m_5 + m_2 m_3 m_4 + m_2 m_3 m_5 + m_2 m_4 m_5 + m_3 m_4 m_5 + m_2 \theta \varepsilon v + m_3 \theta \varepsilon v + m_4 \theta \varepsilon v - g_1 m_1 m_3 - g_1 m_1 m_4 - g_1 m_1 m_5 - g_1 m_3 m_4 - g_1 m_3 m_5 - g_1 m_3 \lambda - g_1 m_4 m_5 - g_1 m_4 \lambda - g_1 m_5 \lambda - g_1 \theta v - \gamma \lambda \theta - g_2 m_1 \sigma - g_2 m_4 \sigma - g_2 m_5 \sigma - g_2 \sigma \lambda
     \)
    \item \(a_4 = m_1 m_2 m_3 m_4 + m_1 m_2 m_3 m_5 + m_1 m_2 m_4 m_5 + m_1 m_3 m_4 m_5 + m_2 m_3 m_4 m_5 + m_2 m_3 \theta \varepsilon v + m_2 m_4 \theta \varepsilon v + m_3 m_4 \theta \varepsilon v - \omega \lambda \theta \gamma_1 - g_1 m_1 m_3 m_4 - g_1 m_1 m_3 m_5 - g_1 m_1 m_4 m_5 - g_1 m_3 m_4 m_5 - g_1 m_3 m_4 \lambda - g_1 m_3 m_5 \lambda - g_1 m_3 \theta \varepsilon v - g_1 m_4 m_5 \lambda - g_1 m_4 \theta \varepsilon v - \gamma m_3 \lambda \theta - \gamma m_4 \lambda \theta - \gamma \sigma \lambda \theta - g_2 m_1 m_4 \sigma - g_2 m_1 m_5 \sigma - g_2 m_4 m_5 \sigma - g_2 m_4 \sigma \lambda - g_2 m_5 \sigma \lambda - g_2 \sigma \theta v\)\\

    \item \(a_5 = m_1 m_2 m_3 m_4 m_5 + m_2 m_3 m_4 \theta v - \phi \sigma \lambda \theta \gamma_1 - \gamma m_3 m_4 \lambda \theta - \gamma \sigma \lambda \theta - m_3 \omega \lambda \theta \gamma_1 - g_1 m_1 m_3 m_4 m_5 - g_1 m_3 m_4 m_5 \lambda - g_2 m_1 m_4 m_5 \sigma - g_2 m_4 m_5 \sigma \lambda - g_1 m_3 m_4 \theta v - g_2 m_4 \sigma \theta v
     \)\\
    
\end{itemize}
Based on the Routh-Hurwitz criteria, the real part of Equation (4.20) is negative if \(a_1 > 0, a_5 > 0, a_1 a_2 - a_3 > 0, a_3 (a_1 a_2 - a_3) + a_1 (a_5 - a_1 a_4) > 0\), and \((a_1 a_2 - a_3)(a_3 a_4 - a_2 a_5) - (a_5 - a_1 a_4)^2 > 0\).

\subsection{Numerical Simulation}


In this section, the simulation of Meningitis spread is conducted with consideration of vaccination, screening, and treatment factors. The model is solved using the fourth-order Runge-Kutta method with a step size of $h=0.1$, implemented in Matlab. The parameters used in the model are listed in Table 1. The simulation is divided into several stages. First, Model (3.8) is simulated using the parameters in Table 3. In the second stage, the simulation is performed by adjusting the parameter value to $\beta=4.88$. This simulation aims to evaluate the effect of infected individuals on the spread of Meningitis within the population. Subsequently, the simulation focuses on intervention strategies to understand the impact of interventions on the spread of Meningitis.

\begin{table}[h!]
\centering
\begin{tabular}{|c|c|c|}
\hline
Parameter & Description & Value \\
\hline
$\mu$ & Natural birth/death rate & 0.02 \\
$\beta$ & Effective contact rate & 5.88 \\
$\eta$ & Per capita infection rate by carriers & [0-1] \\
$\eta_1$ & Per capita infection rate by ill individuals & [0-1] \\
$v$ & Vaccination rate against Meningitis & [0-1] \\
$\theta$ & Loss of immunity & 0.04-2 \\
$\omega$ & Tracing rate & [2-3] \\
$\sigma$ & Rate of progression from C to I & 0.1-0.52 \\
$\phi$ & Treatment rate & 0.3 \\
$\delta$ & Disease-induced mortality & 0.05-0.5 \\
$\gamma$ & Natural recovery rate from disease & 0.1-0.15 \\
$\gamma_1$ & Recovery rate with treatment & 0.05-0.2 \\
\hline
\end{tabular}
\caption{Parameters in Simulation}
\label{tabel2}
\end{table}


\begin{itemize}
    \item {Endemic Equilibrium Simulation}
\end{itemize}

\begin{figure}[htbp]
        \centering
        \includegraphics[width=0.75\linewidth]{EE.png}
        \caption{Solution of Model (3.8) using Parameters in \ref{tabel2}}
        \label{fig:sim_EE}
    \end{figure}
       

Figure 4 shows the simulation results using the parameters in Table \ref{tabel2}. The simulation yields \( R_0 = 1.1228 > 1 \). The value of \( R_0 \) indicates that Meningitis is still spreading within the population. The disease-free and endemic equilibrium points obtained are \( E_0 = (0.55002, 0, 0, 0, 0.44998) \) and \( E_1 = (0.48987, 0.00267, 0.00094, 0.01040, 0.46777) \), respectively. In Figure 4, the simulation shows that the model is locally stable at the endemic point \( E_1 \).

\begin{figure}[htbp]
    \centering
        \includegraphics[width=0.75\linewidth]{DFE.png}
        \caption{Solution of Model (8) using Parameters in Table 1 except \( \beta = 4.88 \)}
        \label{fig:sim\( \beta = 4.88 \) }
    \end{figure}

Figure 5 presents the numerical simulation results with the parameter \( \beta \) adjusted to \( \beta = 4.88 \), indicating a reduction in the interaction rate between infected and susceptible individuals. This adjustment results in a decreased \( R_0 \) value of \( R_0 = 0.9318 < 1 \). According to the simulation results, the model is locally stable at the equilibrium point \( E_0 = (0.55002, 0, 0, 0, 0.44998) \).

\begin{itemize}
    \item {The Impact of Vaccination}
\end{itemize}

\begin{figure}[htbp]
    \centering
    
        \includegraphics[width=0.75\linewidth]{Vaccination Proportion.png}

    \caption{Solution of Model (8) using Parameters in Table 1 except \( v = 0.2 \)}
\end{figure}

In this section, numerical simulation is used to examine the effect of vaccination proportion on the spread of Meningitis. The initial simulation parameter values result in \( R_0 = 1.1228 \). The vaccination proportion parameter in Table \ref{tabel2} is changed to 0.2. The simulation shows that increasing the vaccination proportion decreases the \( R_0 \) value to \( R_0 = 0.7024 \) with a disease-free equilibrium point \( E_0 = (0.37934, 0, 0, 0, 0.62066) \). The effect of vaccination proportion on the spread of Meningitis is shown in Figure 6.

\begin{itemize}
    \item {The Impact of Screening}
\end{itemize}

The numerical simulation in this section aims to show the effect of early screening on individuals carrying the Meningitis bacteria (carrier compartment). The initial simulation uses the parameter value \( \omega = 2.5 \), resulting in \( R_0 = 1.1228 \). Figure 7 shows the numerical solution results with the screening parameter value changed to \( \omega = 3 \). This change means an acceleration in the screening rate for individuals. This adjustment results in a decreased \( R_0 \) value of \( R_0 = 0.9586 < 1 \) with a disease-free equilibrium point \( E_0 = (0.55002, 0, 0, 0, 0.44998) \).

\begin{figure}[htbp]
    \centering

        \includegraphics[width=0.75\linewidth]{Screening Proportion.png} 
        \caption{Solution of Model (8) using Parameters in Table 1 except \( \omega = 3 \)}
\end{figure}

\section{Conclusion}

This study focuses on the spread of Meningitis in a population, considering two mitigation factors: vaccination and screening. The constructed model is divided into five compartments according to the clinical classification of individuals. Two equilibrium points are obtained in this model: the disease-free equilibrium point and the endemic equilibrium point. Using the next-generation matrix method, the basic reproduction number \( R_0 \) is obtained. If \( R_0 < 1 \), the disease-free equilibrium point \( E_0 \) is locally stable. According to the Routh-Hurwitz criteria, the endemic equilibrium point is locally asymptotically stable if \( R_0 > 1 \) and meets the Routh-Hurwitz criteria. Sensitivity analysis shows that the parameters \( v \) and \( \omega \) significantly impact the reduction of \( R_0 \) and the stability of the equilibrium points. This is verified through numerical simulation using the fourth-order Runge-Kutta method with Matlab. The simulation results show that increasing the values of parameters \( v \) and \( \omega \), which means increasing the vaccination proportion and the screening rate, will impact the transmission of Meningitis in the population.

\section{Funding}

This research is funded by the internal Research grant \textit{Penelitian Dosen Pemula (PDP)} 2024, organized by the LPPM Institut Teknologi Bacharuddin Jusuf Habibie.



\begin{thebibliography}{99}

\bibitem{ref1}
D. Aldila \textit{et al.}, ``A mathematical study on the spread of COVID-19 considering social distancing and rapid assessment: The case of Jakarta, Indonesia,'' \textit{Chaos, Solitons \& Fractals}, vol. 139, p. 110042, Oct. 2020, doi: 10.1016/j.chaos.2020.110042.

\bibitem{ref2}
Attaullah and M. Sohaib, ``Mathematical modeling and numerical simulation of HIV infection model,'' \textit{Results in Applied Mathematics}, vol. 7, p. 100118, Aug. 2020, doi: 10.1016/j.rinam.2020.100118.

\bibitem{ref3}
A. Abadi, M. Fakhruddin, R. Artiono, and B. P. Prawoto, ``Measles Transmission Model with Vaccination and Hospitalization Treatments,'' \textit{CBMS}, vol. 3, no. 2, pp. 127–134, May 2021, doi: 10.5614/cbms.2020.3.2.4.

\bibitem{ref4}
F. K. Alalhareth, U. Atta, A. H. Ali, A. Ahmad, and M. H. Alharbi, ``Analysis of Leptospirosis transmission dynamics with environmental effects and bifurcation using fractional-order derivative,'' \textit{Alexandria Engineering Journal}, vol. 80, pp. 372–382, Oct. 2023, doi: 10.1016/j.aej.2023.08.063.

\bibitem{ref5}
S. Adeyemo, A. Sangotola, and O. Korosteleva, ``Modeling Transmission Dynamics of Tuberculosis–HIV Co-Infection in South Africa,'' \textit{Epidemiologia}, vol. 4, no. 4, pp. 408–419, Oct. 2023, doi: 10.3390/epidemiologia4040036.

\bibitem{ref6}
M. Fakhruddin \textit{et al.}, ``Investigation of a measles transmission with vaccination: a case study in Jakarta, Indonesia,'' \textit{MBE}, vol. 17, no. 4, pp. 2998–3018, 2020, doi: 10.3934/mbe.2020170.

\bibitem{ref7}
M. Kizito and J. Tumwiine, ``A Mathematical Model of Treatment and Vaccination Interventions of Pneumococcal Pneumonia Infection Dynamics,'' \textit{Journal of Applied Mathematics}, vol. 2018, pp. 1–16, 2018, doi: 10.1155/2018/2539465.

\bibitem{ref8}
M. R. Nisardi, K. Kasbawati, K. Khaeruddin, A. Robinet, and K. Chetehouna, ``Fractional Mathematical Model of Covid-19 with Quarantine,'' \textit{InPrime:Ind. Jour. Pure. Applied Math}, vol. 4, no. 1, pp. 33–48, Apr. 2022, doi: 10.15408/inprime.v4i1.23719.

\bibitem{ref9}
P. Widyastuti, H. N. Utami, and M. F. Anugrah, ``Meningitis Bakterial: Epidemiologi, Patofisiologi, Penatalaksanaan,'' vol. 2, no. 2, 2023.

\bibitem{ref10}
I. Abdullahi Baba \textit{et al.}, ``Analysis of meningitis model: A case study of northern Nigeria,'' \textit{AIMS Bioengineering}, vol. 7, no. 4, pp. 179–193, 2020, doi: 10.3934/bioeng.2020016.

\bibitem{ref11}
S. Sulma, S. Toaha, and K. Kasbawati, ``Stability Analysis of Mathematical Models of the Dynamics of Spread of Meningitis with the Effects of Vaccination, Campaigns and Treatment,'' \textit{Jurnal Matematika, Statistika dan Komputasi}, vol. 17, pp. 71–81, Aug. 2020, doi: 10.20956/jmsk.v17i1.10031.

\bibitem{ref12}
C. Türkün, M. Gölgeli, and F. M. Atay, ``A mathematical interpretation for outbreaks of bacterial meningitis under the effect of time-dependent transmission parameters,'' \textit{Nonlinear Dyn}, vol. 111, no. 15, pp. 14467–14484, Aug. 2023, doi: 10.1007/s11071-023-08577-6.

\bibitem{ref13}
S. S. Musa, S. Zhao, N. Hussaini, A. G. Habib, and D. He, ``Mathematical modeling and analysis of meningococcal meningitis transmission dynamics,'' \textit{Int. J. Biomath.}, vol. 13, no. 01, p. 2050006, Jan. 2020, doi: 10.1142/S1793524520500060.

\bibitem{ref14}
M. J. F. Martínez, E. G. Merino, E. G. Sánchez, J. E. G. Sánchez, A. M. D. Rey, and G. R. Sánchez, ``A Mathematical Model to Study the Meningococcal Meningitis,'' \textit{Procedia Computer Science}, vol. 18, pp. 2492–2495, 2013, doi: 10.1016/j.procs.2013.05.426.

\bibitem{ref15}
H. E. Bashir, M. Laundy, and R. Booy, ``Diagnosis and treatment of bacterial meningitis,'' 2003.

\bibitem{ref16}
Direktorat Jenderal Pencegahan dan Pengendalian Penyakit Kemenkes RI, ``FAQ MENINGITIS MENINGOKOKUS,'' 2023. Accessed: Feb. 10, 2024. [Online]. Available: \url{https://infeksiemerging.kemkes.go.id/document/frequently-asked-questions-meningitis-meningokokus/}

\bibitem{ref17}
J. K. K. Asamoah, F. Nyabadza, B. Seidu, M. Chand, and H. Dutta, ``Mathematical Modelling of Bacterial Meningitis Transmission Dynamics with Control Measures,'' \textit{Computational and Mathematical Methods in Medicine}, vol. 2018, pp. 1–21, 2018, doi: 10.1155/2018/2657461.

\bibitem{ref18}
K. G. Varshney and Y. K. Dwivedi, ``Mathematical Modelling of Influenza-Meningitis under the Quarantine effect of influenza,'' 2021.

\bibitem{ref19}
B. S. Kotola and T. T. Mekonnen, ``Mathematical model analysis and numerical simulation for codynamics of meningitis and pneumonia infection with intervention,'' \textit{Sci Rep}, vol. 12, no. 1, p. 2639, Feb. 2022, doi: 10.1038/s41598-022-06253-0.

\bibitem{ref20}
F. Keshtkar, G. H. Erjaee, and M. Boutefnouchet, ``On stability of equilibrium points in nonlinear fractional differential equations and fractional Hamiltonian systems,'' \textit{Complexity}, vol. 21, Jul. 2014, doi: 10.1002/cplx.21581.

\bibitem{ref21}
N. Chitnis, J. Hyman, and J. Cushing, ``Determining Important Parameters in the Spread of Malaria Through the Sensitivity Analysis of a Mathematical Model,'' \textit{Bulletin of mathematical biology}, vol. 70, pp. 1272–1296, Jan. 2008, doi: 10.1007/s11538-008-9299-0.

\end{thebibliography}
\end{document}
