Media campaigns,early diagnosis,isolation and treatment on bacterial meningitis outbreak prevention: A modeling study

Abstract

Bacterial meningitis is a severe infection affecting the protective membranes of the brain and spinal cord, with rapidly worsening symptoms that can lead to life-threatening complications. This study presents an autonomous deterministic epidemic model, $S{I}_a{I}_{s}MRS$ , to explore the dynamics of bacterial meningitis in a community implementing control strategies like media coverage, early diagnosis, isolation, and treatment. We adjust the transmission probability based on media coverage and calculate the effective reproduction number, ${\mathcal{R}}_{ef}$ , which includes contributions from asymptomatic ${\mathcal{R}}_{ef}^a$ and symptomatic individuals ${\mathcal{R}}_{ef}^s$ . We derive the basic reproduction number, ${\mathcal{R}}_0$ , to characterize initial infection spread and analyze the local stability of infection-free and endemic equilibria, using the Routh-Hurwitz criteria. For global stability, we apply Castillo-Chavez method for the infection-free equilibrium and the Lyapunov functional technique for the endemic equilibrium, after a uniform persistence study. A local sensitivity analysis evaluates the impact of each parameter on the threshold dynamical parameters ${\mathcal{R}}_{ef}^a$ and ${\mathcal{R}}_{ef}^s$ . We also explore an optimal control problem using Pontryagin's maximum principle. The paper concludes with numerical simulations that bridge theoretical and numerical findings.

Introduction

The exploration of bacterial meningitis has gained notable attention in the field of mathematical epidemiology, driven by its severity and the rapid transmission dynamics within populations.

Meningitis, also known as cerebrospinal meningitis, is a serious infection marked by the inflammation of the meninges (the protective membranes surrounding the brain and spinal cord). Key symptoms include heightened sensitivity to light, intense headaches, neck stiffness, difficulty in feeding, elevated fever, lethargy, persistent crying in infants, and nausea. The impact of meningitis remains a major public health challenge, even though there have been a considerable effort in reducing the spread, through vaccination and treatment. Meningitis and sepsis are estimated to cause more deaths in children under the age of five than malaria. ¹ The annual global disease induced deaths due to meningitis is about 380,000.^1,2 Recovered individuals can present severe abnormalities with appreciable social and economic cost, these sequelae may occur in 10% to 20% of survivors.^1,2 The types of infectious meningitis are: Fungal, Parasitic, Amebic, Bacterial, and Viral. The most severe forms are bacterial and viral meningitis. Viral meningitis is serious but mostly self-limiting, whereas bacterial meningitis usually requires prompt medical attention. ²

Bacterial meningitis can be transmitted from one individual to another, although certain pathogens, such as Listeria monocytogenes, may also spread through contaminated food sources. If not treated promptly, this infection can be fatal, with mortality rates ranging from 5% to 10%, often occurring within one to two days of symptoms onset. Without appropriate treatment, the death rate can reach up to 50%. Between 1991 and 2010, nearly one million suspected cases were reported, predominantly in countries within the African Meningitis Belt, resulting in approximately 100, 000 deaths, of which 79, 292 of these cases occurred in 2009 with 4, 288 deaths in 14 African countries under the enhanced surveillance.^2,3 During the 2014 and 2013 epidemiological season there were 14, 317 suspected cases with, 1, 304 deaths in 19 African countries, 12, 464 cases with 1, 131 deaths in 18 African countries, this is relatively considered as a mild meningitis epidemic period, since 2005 (see^1–3 for a brief report). The most common forms of bacterial meningitis are meningococcal meningitis and Haemophilus influenzae type b (Hib). Hib is responsible for large outbreaks in sub-Saharan Africa and primarily affects people between the ages of 1 and 30..^1,4 An outbreak of meningococcal meningitis resulted in 11, 000 deaths in the slums of Brazil, during the period of 1974−1978, and left about 75, 000 persons with permanent neurological complications.^4,5 In the Sahel region of West Africa, an outbreak resulted in 25, 000 deaths during the 1996−1997 epidemic season.^1,6 Also, there are about 1400–3000 annually reported cases of invasive meningococcal disease in U.S..^2,4,7 Meningococcal meningitis is caused by the bacterium Neisseria meningitidis. Twelve types of Neisseria meningitidis have been identified and called serogroups. Of these 12 serogroups, six (A, B, C, W, X, and Y) can cause an epidemic. ¹ There are some vaccines available against the three main bacterial meningitis (Nerisseria meningitidis, Haemophilus influenzae type b and Streptococcus pneumoniae). ^{1 8–10} However, there is no vaccine yet against group B streptococcus (Streptococcus agalactiae).^1,4,9,10 In May 2018, the WHO launched the “Defeating Meningitis by 2030” initiative, outlining a strategic framework to combat meningitis outbreaks. This framework focuses on prevention and epidemic control, improved diagnosis and treatment including the development of rapid diagnostic tests-disease surveillance, advocacy and information dissemination, as well as support and aftercare for survivors.^1,11,12 Also, studies on meningococcal disease control, including the WHO's 2019 guidelines on meningitis outbreaks (reporting 60−85% transmission reduction via droplet precautions and 24-h isolation post-antibiotics) and Borrow et al.'s 2017 retrospective analysis (showing 75% efficacy in preventing UK secondary cases through chemoprophylaxis and isolation compliance), highlight the strong impact of these interventions. ¹³ Furthermore, media coverage plays a vital role in the fight against epidemic or pandemic of communicable diseases, especially in the earliest moments of the outbreak.^11,12,14,15 It plays an educational role helping to raise the level of awareness by broadcasting information about the epidemic spreading mechanism and how to behave properly. Therefore, various researchers have developed mathematical models to assess how media coverage influences the transmission dynamics of infectious diseases, utilizing different media functions based on available data. In the context of modeling the impact of media coverage on disease transmission dynamics, two primary approaches are prevalent.^16–18 One method involves treating the transmission term or effective contact rate as a variable function rather than a constant, allowing feedback to adjust transmission likelihood based on the number of infected individuals. For example, Liu et al. ¹² utilized a complex network model in their EIH framework, where H represents hospitalized individuals, and defined the transmission coefficient as $\beta (E,H,I)=e^{-{\alpha }_{1}E-{\alpha }_{2}I-{\alpha }_{3}H}$ . Additionally, G. P. Sahu et al. ¹⁹ explored the dynamics of an SEQIHRS epidemic model incorporating media coverage, quarantine, and isolation measures within a community. In their proposed model, they have used a media-induced transmission rate of the form $\tilde{\beta }{e}^{-m_1\frac{\tilde{I}}{\tilde{N}}-m_2\frac{\tilde{H}}{\tilde{N}}}$ , where $\tilde{\beta }$ is the contact rate between susceptible and infectious persons and the coefficients $m_1$ and $m_2$ represent the non-pharmaceutical interventions of media coverage targeting infectious $\tilde{I}$ and isolated $\tilde{H}$ individuals respectively. Furthermore, P. Song and Y. Xiao ¹¹ investigated the impact of media on an epidemic model by incorporating the following media function $\beta (I)={\beta }_1-\frac{{\beta }_{2}I}{m+I}$ , I representing infected individuals and m a positive constant. The other approach of considering the feedback effect from social media is to treat it as an independent state variable within the model as $f(M)=e^{-\mu M}$ in Ref. ²⁰ , where M denotes the magnitude of media coverage, $f(M)=\frac{M}{m+M}$ in Ref. ²¹ or $f(A)=r\left(1-\theta \frac{A}{\omega +A}\right)I$ proposed by AK Misra et al., ¹⁵ A denoting the awareness of the population. In light of these previous works, we propose an optimal control model for bacterial meningitis which is treatable to study the impact of the media coverage, early diagnosis, isolation and treatment in both the carrier and infected stages of the disease herein represented by $u_1$ and $u_2$ , respectively, and the impact of the disease dynamics within hospitals. We utilize the media function defined as $\beta (I_s)={\beta }_s{e}^{-m{I}_s}{I}_s$ and in this expression, ${\beta }_s$ denotes the effective contact rate between susceptible S and symptomatic infectious $I_s$ , while the non-negative parameter m quantifies the psychological effects of media coverage on symptomatic individuals $I_s$ . This formulation captures how media reporting can shape public perception and behavior in response to disease prevalence.

As part of efforts to combat bacterial meningitis, Asamoah et al. ⁴ proposed an SCIRS model of bacterial meningitis transmission incorporating a standard incidence function and a nonlinear recovery rate. Their study focused on the impact of antibiotic efficacy and vaccine use. They further conducted an optimal control analysis under endemic conditions in resource-limited settings characterized by a limited supply of antibiotics and hospital beds by investigating the effects of inoculation. Similarly, Peter et al. ²² developed a mathematical model of the SVEIQRS type using ordinary differential equations (ODEs), which they extended to fractional order to explore the dynamics of meningitis. Their work examined how vaccination and treatment influence transmission dynamics as functions of various fractional parameter values.

Belay et al. ²³ introduced a deterministic $SVCI{D}_{r}RS$ -type model to perform an optimal control study of bacterial meningitis dynamics, focusing on the cost-effectiveness of three control measures: prevention, treatment, and screening. They evaluated different combined control strategies to assess their impact on disease spread and associated costs. Simulation results indicated that maximizing prevention and screening efforts constitutes the most effective strategy. In addition, Sahar et al. ²⁴ addressed the optimal control of meningococcal meningitis transmission using an $SCIRS$ -type model with a linear incidence function. They employed two control measures (vaccination and treatment) to minimize infections, and performed numerical simulations to demonstrate the implementation of these strategies.

Moreover, Amar et al. ²⁵ studied a deterministic SCITRS-type model of meningitis transmission, considering vaccination and screening as intervention methods. Using the fourth-order Runge-Kutta method, they confirmed the steady-state stability results and showed that vaccination and screening significantly affect meningitis transmission dynamics. Other mathematical studies have also been conducted within this framework, addressing various objectives (e.g. Refs.^21,26).

Using the first approach of integrating media coverage so that the force of infection induced by asymptomatic carriers is given by ${\beta }_s{e}^{-m{I}_s}{I}_s$ and taking into account most of the mechanisms involved in the current dynamics of meningitis, this study introduces a novel compartmental structure to advance meningitis transmission modeling. This allows us to assess the impact of large-scale awareness-raising efforts through media coverage, early diagnosis, isolation of asymptomatic carriers and symptomatic infectious individuals for treatment, self-medication, and recovery through hospital treatment on the dynamics of the disease. Our optimal control model takes into account early diagnosis and treatment administered after isolation in order to highlight the impact of good monitoring of the disease dynamics by health structures, in order to limit its spread. Furthermore, as mathematical modeling can provide insight into certain measures using qualitative analyses and numerical simulations, we hereby aim to evaluate the effectiveness of WHO-recommended interventions^1,8 in controlling meningitis spread and meet the target of eradication by 2030. Hence, we believe that this research will give some further insight on how to effectively control the disease, using the media, early diagnosis and treatment.

Model development

In this study, we consider the disease dynamics in a homogeneous population that is structured in five compartments: Susceptible $S(t)$ , Asymptomatic carriers $I_a(t)$ , Symptomatic infectious individuals $I_s(t)$ , group of individuals receiving medical care in isolation denoted by $M(t)$ and Recovered (self-medication or hospital based treatment) $R(t)$ individuals. Therefore, $N(t)=S(t)+I_a(t)+I_s(t)+M(t)+R(t)$ represents the total population at any time t. We assume that, early diagnosed carrier and ill individuals enter into the hospitalized group at the rate $u_1$ and $\gamma$ is the effective diagnosis rate at the asymptomatic stage, ${\gamma }_1$ is the treatment efficacy. Individuals in the isolated for treatment (self-medication or hospital-based prescriptions) are treated at the rate $u_2$ and are assumed to move to the recovery stage, $\omega$ is the rate of natural recovery at the carrier state, and $\alpha$ is the rate at which carrier individual progress to the ill compartment. The assumed rate of natural death in all the compartments is denoted by $\mu$ and, ${\delta }_s$ and ${\delta }_h$ are the disease induced death rate at the symptomatic infectious and hospitalized compartments, respectively. It is known that recovered individuals do not attain permanent immunity against bacterial meningitis after recovery from infection,^22–25,27 hence we assume $\vartheta$ as the rate of regaining susceptibility after recovery. The influx of individuals into the population is assumed to join the susceptible compartment at a rate $\mathrm{\Lambda }$ , which represents the recruitment rate that includes both new births and immigrants. We also assume that the hospitalized individuals are perfectly isolated so they have no chance to infect people. The parameters ${\beta }_a$ and ${\beta }_s$ represent the effective contact rate of susceptible population with asymptotic carriers and symptomatic infected people, respectively. However, the coefficient m represents the impact of media coverage to the contact for symptomatic infectious individuals. The proposed model is illustrated by the schematic flow chart presented in Figure 1.

Figure 1.

Schematic flow chart of the proposed model for bacterial meningitis with the incidence function given by $f(I_a,I_s)={\beta }_s{e}^{-m{I}_s}{I}_s+{\beta }_a{I}_a$ .

Using the aforementioned descriptions and assumptions, and by applying the input–output balance from the flow chart of Figure 1, the dynamics of the disease are described by the following system of ODEs:

$$\{\begin{array}{c}{\displaystyle \frac{dS}{dt}=\mathrm{\Lambda }+\vartheta R-{\beta }_s{e}^{-m{I}_s}S{I}_s-{\beta }_{a}S{I}_a-\mu S}\\ {\displaystyle \frac{d{I}_a}{dt}={\beta }_s{e}^{-m{I}_s}S{I}_s+{\beta }_{a}S{I}_a-(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )I_a,}\\ {\displaystyle \frac{d{I}_s}{dt}=(1-\gamma u_1)\alpha I_a-({\delta }_s+u_1+\mu )I_s,}\\ {\displaystyle \frac{dM}{dt}=\gamma u_1{I}_a+u_1{I}_s-({\delta }_h+{\gamma }_1{u}_2+\mu )M,}\\ {\displaystyle \frac{dR}{dt}={\gamma }_1{u}_{2}M+\omega I_a-(\vartheta +\mu )R,}\end{array}$$

(1)

with initial conditions

$$\begin{array}{rl}S(0)=& S^0>0,I_a(0)=I_a^0\ge 0,I_s(0)=I_s^0\ge 0,M(0)=M^0\ge 0,R(0)=R^0\ge 0,\\ & N(0)=N^0\ge 0.\end{array}$$

Descriptions of all the parameters of the model are subscribed in Table 1.

Table 1.

Model parameter description.

Parameter	Description	Unit
$\mathrm{\Lambda }$	Recruitment rate (births and immigrants)	humans/day
${\beta }_s$	Rate of effective contact between susceptible and symptomatic ( $I_s$ )	day−1
${\beta }_a$	Rate of effective contact between susceptible and asymptomatic ( $I_a$ )	day−1
$\mu$	Rate of natural mortality	day−1
${\delta }_s$	Disease death rate of symptomatic infectious ( $I_s$ )	day−1
${\delta }_h$	Disease death rate among hospitalized patients ( $M$ )	day−1
$u_1$	Early diagnosis control parameter	day−1
$\gamma$	Early diagnosis efficacy parameter of asymptomatic carriers	day−1
$u_2$	Treatment rate of hospitalized	day−1
${\gamma }_1$	Treatment efficacy	day−1
$\omega$	Natural recovery rate of symptomatic carriers	day−1
$\alpha$	Transition rate from the compartment ( $I_a$ ) to compartment ( $I_s$ )	day−1
$\vartheta$	Immunity loss rate	day−1
$1/\vartheta$	Average waning period of disease-induced immunity	day
$m$	Influence of media coverage on symptomatic infection contact	day−1

Mathematical analysis of the model

Results on positivity invariance and boundedness

Note that

$${\displaystyle \frac{dN}{dt}=\mathrm{\Lambda }-\mu N-{\delta }_s{I}_s-{\delta }_{h}M}$$

(2)

This study focuses solely on solutions with initial conditions within the following biologically viable set, ensuring the existence, uniqueness, and continuation of solutions:

$$\mathrm{\Gamma }=\left\{(S,I_a,I_s,M,R)\in {\mathbb{R}}_{+}^5:0\le S,I_a,I_s,M,R,N\le {\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\right\}$$

Furthermore, it is easy to show the boundedness and the positively invariant of the region Γ with respect to system (1).

Proposition 3.1. All solution trajectories of system

(1) that begin within the set $\mathrm{\Gamma }$ will either approach, enter, or remain within the interior of $\mathrm{\Gamma }$ .

Proof.

Let ${\mathbb{R}}_{+}^5=\{(S,I_a,I_s,M,R)\in {\mathbb{R}}^5:S\ge 0,I_a\ge 0,I_s\ge 0,M\ge 0,R\ge 0\}$ denotes the non-negative cone in the Euclidean space of dimension 5. From system (1),

$$\begin{array}{rl}{\displaystyle \frac{dS}{dt}{|}_{S=0}=}& \mathrm{\Lambda }+\vartheta R>0,{\displaystyle \frac{d{I}_a}{dt}{|}_{I_a=0}={\beta }_s{e}^{-m{I}_s}S{I}_s>0,{\displaystyle \frac{d{I}_s}{dt}{|}_{I_s=0}=(1-\gamma u_1)\alpha I_a>0,}}\\ {\displaystyle \frac{dM}{dt}{|}_{M=0}=}& \gamma u_1{I}_a+u_1{I}_s>0,{\displaystyle \frac{dR}{dt}{|}_{R=0}={\gamma }_1{u}_{2}M+\omega I_a>0,}\end{array}$$

and $S(t),I_a(t),I_s(t),M(t),R(t)$ are continuous functions of t. Therefore, the vector field on each bounding hyperplane of ${\mathbb{R}}_{+}^5$ directs inward, ensuring that all solution trajectories starting in ${\mathbb{R}}_{+}^5$ will remain within it at all time. As a result, ${\mathbb{R}}_{+}^5$ is positively invariant for system (1). Moreover, from equation (2), we have $\frac{dN}{dt}=\mathrm{\Lambda }-\mu N-{\delta }_s{I}_s-{\delta }_{h}M$ and then have $\frac{dN}{dt}<\mathrm{\Lambda }-\mu N.$ By applying Birkhoff's and Rota's theorem on differential inequalities, ¹⁵ we find that as $t\to \mathrm{\infty },$ the inequality $0\le N(t)\le \frac{\mathrm{\Lambda }}{\mu }$ holds. This indicates that the solution of system (1) is bounded, ensuring that any solution originating from $\mathrm{\Gamma }$ will remain within $\mathrm{\Gamma }$ .^16–18

Local stability analysis of the disease-free equilibrium

In this section, we calculate the disease-free equilibrium $dfe=\left(\frac{\mathrm{\Lambda }}{\mu },0,0,0,0\right)$ and the effective reproduction number ${\mathcal{R}}_{ef}.$ We also examine the local stability of the $dfe.$ The non-negative matrix $\mathcal{F}$ , representing the new infection terms, and the matrix $\mathcal{V}$ , which encompasses the remaining terms, are defined as follows:

$$\mathcal{F}=\left(\begin{array}{c}{\beta }_s{e}^{-m{I}_s}S{I}_s+{\beta }_{a}S{I}_a\\ \begin{array}{c}0\\ 0\end{array}\end{array}\right)and\mathcal{V}=\left(\begin{array}{c}(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )I_a\\ \begin{array}{c}-(1-\gamma u_1)\alpha I_a+({\delta }_s+u_1+\mu )I_s\\ -\gamma u_1{I}_a-u_1{I}_s+({\delta }_h+{\gamma }_1{u}_2+\mu )M\end{array}\end{array}\right).$$

Linearizing both matrices at the $dfe$ leads to the following matrices:

$$F=\left(\begin{array}{ccc}{\displaystyle \frac{{\beta }_a\mathrm{\Lambda }}{\mu }}& {\displaystyle \frac{{\beta }_s\mathrm{\Lambda }}{\mu }}& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)andV=\left(\begin{array}{ccc}\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu & 0& 0\\ -(1-\gamma u_1)\alpha & {\delta }_s+u_1+\mu & 0\\ -\gamma u_1& -u_1& {\delta }_h+{\gamma }_1{u}_2+\mu \end{array}\right).$$

The effective reproduction number is the largest eigenvalue of the next generation matrix. Therefore,

$${\mathcal{R}}_{ef}={\mathcal{R}}_{ef}^a+{\mathcal{R}}_{ef}^s$$

(3)where

$${\mathcal{R}}_{ef}^a={\displaystyle \frac{{\beta }_a\mathrm{\Lambda }}{\mu (\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )}}$$

and

$${\mathcal{R}}_{ef}^s={\displaystyle \frac{{\beta }_s\mathrm{\Lambda }(1-\gamma u_1)\alpha }{\mu (\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )({\delta }_s+u_1+\mu )}.}$$

From equation (3), it is clear that the threshold quantity for our proposed model system consists of two parts, which means that there are two different transmissions routes; one from the asymptomatic carriers ${\mathcal{R}}_{ef}^a$ and the other from the symptomatic infectious individuals ${\mathcal{R}}_{ef}^s$ . In the realm of infectious disease epidemiology, the concepts of basic and effective reproduction numbers are fundamental to understanding the dynamics of disease transmission. While often used interchangeably, these two metrics have distinct implications for public health policy and control intervention strategies.

The basic reproduction number ${\mathcal{R}}_0$ represents the average number of secondary cases generated by a single infected individual in a completely susceptible population, absent any immunity from past exposures or vaccination. This metric provides a measure of the intrinsic transmissibility of a pathogen, reflecting its ability to spread in a population with no prior immunity.

In contrast, the effective reproduction number ${\mathcal{R}}_{ef}$ represents the average number of secondary cases generated by a single infected individual in a population with some level of immunity, either from past exposures or vaccination, and in the presence of intervention measures including contact tracing, quarantine, or vaccination campaigns. This metric takes into account the impact of these interventions on the transmission dynamics of the disease.

Therefore, setting $\gamma ={\gamma }_1=0$ which means there is no isolation, no early diagnosis, no treatment within the community, one can drop the expression of the basic reproduction number as follows:

$${\mathcal{R}}_0={\mathcal{R}}_0^a+{\mathcal{R}}_0^s$$

(4)where

$${\mathcal{R}}_0^a={\displaystyle \frac{{\beta }_a\mathrm{\Lambda }}{\mu (\omega +\alpha +\mu )}and{\mathcal{R}}_0^s={\displaystyle \frac{{\beta }_s\mathrm{\Lambda }\alpha }{\mu (\omega +\alpha +\mu )({\delta }_s+\mu )}.}}$$

In Figure 2, we show the evolution of ${\mathcal{R}}_0$ and ${\mathcal{R}}_{ef}$ as a function of multivariable.

Figure 2.

(a) Variation of ${\mathcal{R}}_0$ as a function of ${\mathcal{R}}_0^a$ and ${\mathcal{R}}_0^s$ , (b) variation of ${\mathcal{R}}_{ef}$ as a function of ${\mathcal{R}}_{ef}^a$ and ${\mathcal{R}}_{ef}^s$ .

From Theorem 2 in Ref., ²⁸ the following result holds:

Theorem 3.1. The

meningitis-free equilibrium of system (1) exhibits local asymptotic stability when ${\mathcal{R}}_{ef}$ < 1, while it transitions to instability when ${\mathcal{R}}_{ef}$ > 1.

In the sense of epidemiology, this result means that if the effective reproduction number ${\mathcal{R}}_{ef}<1$ , then the influx of a few infected individuals will not generate large outbreaks (and the disease will die out) whereas the disease outbreak will occur if ${\mathcal{R}}_{ef}>1.$

The $dfe$ global stability analysis

To establish the global stability of the disease-free steady state, we use the method developed by Castillo-Chavez et al. ²⁹ For that, we rewrite the system (1) in the following form:

$$\{\begin{array}{c}{\displaystyle \frac{dX}{dt}=\mathcal{K}(X,Y),}\\ {\displaystyle \frac{dY}{dt}={\mathcal{K}}_1(X,Y),\mathcal{G}(X,0)=0,}\end{array}$$

(5)where $X={(S)}^T\in {\mathbb{R}}_{+}$ represents the number of uninfected individuals and $Y=(I_a,I_s,M{)}^T$ $\in {\mathbb{R}}_{+}^3$ denotes the number of infected individuals that comprises asymptomatic carriers, symptomatic infectious and people isolated with medical care. The disease-free equilibrium hereby denoted by $dfe=(X_{0,}0)=\left(\frac{\mathrm{\Lambda }}{\mu },0\right)$ is globally asymptotically stable if the following conditions are satisfied:

(H1) for $\frac{dX}{dt}=\mathcal{K}(X,0),X_0$ is globally asymptotically stable,

(H2) ${\mathcal{K}}_1(X,Y)=\mathcal{M}Y-\mathcal{G}(X,Y),\mathcal{G}(X,Y)\ge 0$ for $(X,Y)\in \mathrm{\Gamma }$ ,

and $\mathcal{M}=D_Y{\mathcal{K}}_1(X_0,0)$ is Metzler stable. Whenever our model system (1) of differential equations describing the disease spreading dynamics within the population obeys hypothesis (H1) and (H2), then the infection-free point $dfe=(X_0,0)$ is a globally asymptotically stable steady state of the mathematical model provided when ${\mathcal{R}}_{ef}<1$ .

This is better expressed by the forthcoming result:

Theorem 3.2. If the effective reproduction number

${\mathcal{R}}_{ef}<1$ , then the disease-free equilibrium point $dfe=(X_0,0)$ is globally asymptotically stable.

Proof.

From system (1), we draw $\mathcal{K}(X,Y)$ and ${\mathcal{K}}_1(X,Y)$ as follows:

$$\mathcal{K}(X,Y)=\mathrm{\Lambda }+\vartheta R-{\beta }_s{e}^{-m{I}_s}S{I}_s-{\beta }_{a}S{I}_a-\mu S$$

(6)and

$${\mathcal{K}}_1(X,Y)=\left(\begin{array}{c}{\beta }_s{e}^{-m{I}_s}S{I}_s+{\beta }_{a}S{I}_a-(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )I_a\\ \begin{array}{c}(1-\gamma u_1)\alpha I_a-({\delta }_s+u_1+\mu )I_s\\ \gamma u_1{I}_a+u_1{I}_s-({\delta }_h+{\gamma }_1{u}_2+\mu )M\end{array}\end{array}\right).$$

(7)

At Y = 0, equation (6) leads to

$$\mathcal{K}(X,0)=\mathrm{\Lambda }-\mu S$$

(8)and using the method of integration factor yields

$$S(t)={\displaystyle \frac{\mathrm{\Lambda }}{\mu }-S(0)e^{-\mu t}.}$$

(9)

Thus, for any value of the initial condition $S(0)$ , it follows that as $t\to \mathrm{\infty },S(t)\to \frac{\mathrm{\Lambda }}{\mu }.$

Therefore, $X_0$ is globally asymptotically stable. Hence, condition (H1) is satisfied.

Now, let us consider the infected compartment of the meningitis given by

$${\mathcal{K}}_1(X,Y)=\mathcal{M}Y-\mathcal{G}(X,Y)$$

(10)where,

$$\mathcal{M}=\left(\begin{array}{ccc}(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )({\mathcal{R}}_{ef}^a-1)& {\displaystyle \frac{{\beta }_s\mathrm{\Lambda }}{\mu }}& 0\\ (1-\gamma u_1)\alpha & -({\delta }_s+u_1+\mu )& 0\\ \gamma u_1& u_1& -({\delta }_h+{\gamma }_1{u}_2+\mu )\end{array}\right)$$

and

$$\mathcal{G}(X,Y)=\left(\begin{array}{c}{\beta }_s(S^0-e^{-m{I}_s}S)I_s+{\beta }_a(S^0-S)I_a\\ \begin{array}{c}0\\ 0\end{array}\end{array}\right).$$

Since ${\mathcal{R}}_{ef}^a<{\mathcal{R}}_{ef}<1$ , it is easy to see that the matrix $\mathcal{M}$ is Metzler stable. Furthermore, for all $(X,Y)\in \mathrm{\Gamma },(S^0-S)>0$ and $(S^0-e^{-m{I}_s}S)>0$ , therefore $\mathcal{G}(X,Y)>0.$

This proves that the condition (H2) is satisfied. Hence, the disease-free equilibrium is globally asymptotically stable whenever ${\mathcal{R}}_{ef}<1$ .

Presence of the endemic steady state and its local stability

Setting the right-hand side of the system (1) to zero, we obtain the endemic equilibrium point as $\mathcal{E}*=(S*,I_a*,I_s*,M*,R*)$ so that

$$\begin{array}{rl}{I}_s*=& {\displaystyle \frac{(1-\gamma u_1)\alpha }{{\delta }_s+u_1+\mu }{I}_a*=\mathfrak{F}(I_a*),S*={\displaystyle \frac{(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )({\delta }_s+u_1+\mu )}{{\beta }_s(1-\gamma u_1)\alpha e^{-m\mathfrak{F}(I_a*)}+{\beta }_a({\delta }_s+u_1+\mu )}}}\\ M*=& {\displaystyle \frac{\gamma u_1({\delta }_s+u_1+\mu )+u_1(1-\gamma u_1)\alpha }{({\delta }_s+u_1+\mu )({\delta }_h+{\gamma }_1{u}_2+\mu )}{I}_a*,}\\ R*=& {\displaystyle \frac{[\gamma u_1({\delta }_s+u_1+\mu )+u_1(1-\gamma u_1)\alpha ]{\gamma }_1{u}_2+\omega ({\delta }_s+u_1+\mu )({\delta }_h+{\gamma }_1{u}_2+\mu )}{({\delta }_s+u_1+\mu )({\delta }_h+{\gamma }_1{u}_2+\mu )(\vartheta +\mu )}{I}_a*,}\end{array}$$

where $I_a*$ satisfies the following equation:

$$\begin{array}{rl}\mathcal{W}(I_a*)=& \mathrm{\Lambda }+{\displaystyle \frac{\vartheta [\gamma u_1({\delta }_s+u_1+\mu )+u_1(1-\gamma u_1)\alpha ]{\gamma }_1{u}_2+\vartheta \omega ({\delta }_s+u_1+\mu )({\delta }_h+{\gamma }_1{u}_2+\mu )}{({\delta }_s+u_1+\mu )({\delta }_h+{\gamma }_1{u}_2+\mu )(\vartheta +\mu )}}\\ & I_a*-{\displaystyle \frac{{\beta }_s\alpha (\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )(1-\gamma u_1)I_a*e^{-m\mathfrak{F}(I_a*)}}{{\beta }_s\alpha (1-\gamma u_1)e^{-m\mathfrak{F}(I_a*)}+{\beta }_a({\delta }_s+u_1+\mu )}}\\ & -{\displaystyle \frac{{\beta }_a(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )({\delta }_s+u_1+\mu )I_a*}{{\beta }_s\alpha (1-\gamma u_1)e^{-m\mathfrak{F}(I_a*)}+{\beta }_a({\delta }_s+u_1+\mu )}}\\ & -{\displaystyle \frac{\mu (\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )({\delta }_s+u_1+\mu )}{{\beta }_s\alpha (1-\gamma u_1)e^{-m\mathfrak{F}(I_a*)}+{\beta }_a({\delta }_s+u_1+\mu )}=0.}\end{array}$$

Indeed,

$$\begin{array}{rl}\mathcal{W}(0)=& \mathrm{\Lambda }+{\displaystyle \frac{\mu (\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )({\delta }_s+u_1+\mu )}{{\beta }_s\alpha (1-\gamma u_1)+{\beta }_a({\delta }_s+u_1+\mu )}}\\ & =\mathrm{\Lambda }\left(1-{\displaystyle \frac{\mu (\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )({\delta }_s+u_1+\mu )}{{\beta }_s\mathrm{\Lambda }\alpha (1-\gamma u_1)+{\beta }_a\mathrm{\Lambda }({\delta }_s+u_1+\mu )}}\right)=\mathrm{\Lambda }\left(1-{\displaystyle \frac{1}{{\mathcal{R}}_{ef}}}\right)>0\end{array}$$

for ${\mathcal{R}}_{ef}>1.$

Besides,

$${\mathcal{W}}^{\mathrm{″}}(I_a*)={\displaystyle \frac{2m{\beta }_s\alpha {(1-\gamma u_1)}^2{e}^{-m\mathfrak{F}(I_a*)}}{\alpha ({\delta }_s+u_1+\mu ){({\beta }_s\alpha (1-\gamma u_1)e^{-m\mathfrak{F}(I_a*)}+{\beta }_a({\delta }_s+u_1+\mu ))}^3}{\mathcal{W}}_1(I_a*)}$$

(11)with

$$\begin{array}{rl}{\mathcal{W}}_1(I_a*)=& \alpha (1-\gamma u_1)(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )e^{-m\mathfrak{F}(I_a*)}\\ & -{\beta }_a({\delta }_s+u_1+\mu )(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )-1.\end{array}$$

Furthermore, we make the following assumption:

(H3) $\mathcal{W}\left(\frac{\mathrm{\Lambda }}{\mu }\right)\le 0and{\mathcal{W}}_1(I_a*)\le 0for{I}_a*\in \left(0,\frac{\mathrm{\Lambda }}{\mu }\right).$

Therefore, $\mathcal{W}$ is a continuous function with $\mathcal{W}(0)>0$ and $\mathcal{W}\left(\frac{\mathrm{\Lambda }}{\mu }\right)\le 0.$ Thus, $\mathcal{W}$ admits a zero within the interval $\left(0,\frac{\mathrm{\Lambda }}{\mu }\right)$ . However, thanks to the assumption (H3), equation (11) gives ${\mathcal{W}}^{\mathrm{″}}(I_a*)\le 0$ which shows the concavity of the function $\mathcal{W}$ on $\left(0,\frac{\mathrm{\Lambda }}{\mu }\right).$ Consequently, we have the uniqueness of the zero of $\mathcal{W}$ on $\left(0,\frac{\mathrm{\Lambda }}{\mu }\right)$ and then the model has a unique endemic equilibrium.

And the following result holds:

Theorem 3.3. If

(H3) is satisfied and if ${\mathcal{R}}_{ef}>1$ , then model (1) has a unique positive endemic equilibrium $\mathcal{E}*=(S*,I_a*,I_s*,M*,R*).$

Theorem 3.4. The endemic equilibrium point

$\mathcal{E}*$ is locally asymptotically stable whenever ${\mathcal{R}}_{ef}>1$ under the condition that $a_0>0,a_4>0$ and the Routh-Hurwitz determinants ${\mathrm{\Delta }}_{1,1},{\mathrm{\Delta }}_{1,2},{\mathrm{\Delta }}_{2,1},{\mathrm{\Delta }}_{2,2},{\mathrm{\Delta }}_{3,1}$ are positive.

Proof.

Let $d_1=\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu ,d_2={\delta }_s+u_1+\mu ,d_3={\delta }_h+{\gamma }_1{u}_2+\mu$ and $d_4=\vartheta +\mu$ . The Jacobian matrix of System (1) at endemic equilibrium $\mathcal{I}(\mathcal{E}*)$ is given by

$$\mathcal{I}(\mathcal{E}*)=(\begin{array}{ccccc}-{\beta }_s{e}^{-m{I}_s*}{I}_s*-{\beta }_a{I}_a*-\mu & -{\beta }_{a}S*& m{\beta }_s{e}^{-m{I}_s*}S*I_s*-{\beta }_s{e}^{-m{I}_s*}S*& 0& \vartheta \\ {\beta }_s{e}^{-m{I}_s*}{I}_s*+{\beta }_a{I}_a*& {\beta }_{a}S*-d_1& -m{\beta }_s{e}^{-m{I}_s*}S*I_s*+{\beta }_s{e}^{-m{I}_s*}S*& 0& 0\\ 0& \alpha (1-\gamma u_1)& -d_2& 0& 0\\ 0& \gamma u_1& u_1& -d_3& 0\\ 0& \omega & 0& {\gamma }_1{u}_2& -d_4\end{array})$$

and its eigenvalues are the roots of the following polynomial equation:

$${\chi }^5+a_4{\chi }^4+a_3{\chi }^3+a_2{\chi }^2+a_1\chi +a_0=0,$$

(12)where

$$\begin{array}{rl}{a}_4=& {\beta }_s{e}^{-m{I}_s*}{I}_s*+{\beta }_a{I}_a*+\mu +\sum _{i=1}^4{d}_i,\\ a_3=& {\beta }_s{e}^{-m{I}_s*}\sum _{i=1}^4{d}_i+\alpha (1-\gamma u_1)(m{I}_s*-1){\beta }_s{e}^{-m{I}_s*}S*+(\mu +{\beta }_a{I}_a*)\left(\sum _{i=3}^4{d}_i-\sum _{i=1}^2{d}_i\right)\\ & +{\beta }_a^{2}S*(I_a*+I_s*)-{\beta }_{a}S*\sum _{i=2}^4{d}_i+\mu {\beta }_{a}S*+d_1{d}_2+d_1{d}_3+d_1{d}_4+d_2{d}_3+d_2{d}_4+d_3{d}_4,\\ a_2=& ({\beta }_s{e}^{-m{I}_s*}{I}_s*+{\beta }_a{I}_a*)\left[{\beta }_{a}S*\left(1+\sum _{i=3}^4{d}_i\right)+\prod _{i=1}^2{d}_i-\sum _{i=1}^2{d}_i-\vartheta \omega +1\right]+\mu \\ & +\alpha (1-\gamma u_1)(m{I}_s*-1){\beta }_s{e}^{-m{I}_s*}S*-{\beta }_{a}S*(1+\mu )(1+d_2)+(1+\mu )\left(\sum _{i=1}^2{d}_i+\prod _{i=1}^2{d}_i\right)\\ a_1=& \sum _{i=3}^4{d}_i\left[\prod _{i=1}^2{d}_i({\beta }_s{e}^{-m{I}_s*}{I}_s*+{\beta }_a{I}_a*+\mu \alpha (1-\gamma u_1)(m{I}_s*-1){\beta }_s{e}^{-m{I}_s*}S*)-\mu d_2{\beta }_{a}S*+\mu \prod _{i=1}^2{d}_i\right]\\ & +\prod _{i=3}^4{d}_i\left[({\beta }_s{e}^{-m{I}_s*}{I}_s*+{\beta }_a{I}_a*+\mu )\sum _{i=1}^2{d}_i+\alpha (1-\gamma u_1)(m{I}_s*-1){\beta }_s{e}^{-m{I}_s*}S*-{\beta }_{a}S*(\mu +d_2)+\prod _{i=1}^2{d}_i\right]\\ & -({\beta }_s{e}^{-m{I}_s*}{I}_s*+{\beta }_a{I}_a*)\left(\vartheta \omega \sum _{i=1}^2{d}_i+\gamma {\gamma }_1{u}_1{u}_2\vartheta \right)\\ a_0=& ({\beta }_s{e}^{-m{I}_s*}{I}_s*+{\beta }_a{I}_a*)\left(d_1\prod _{i=3}^4{d}_i-\vartheta \omega \prod _{i=2}^3{d}_i-\vartheta {\gamma }_1{u}_1{u}_2\alpha (1-\gamma u_1)-d_2\vartheta \gamma {\gamma }_1{u}_1{u}_2\right)\end{array}$$

Now, applying Routh-Hurwitz criterion one can see that the following conditions are necessary for the existence of negative roots or negative real parts of the roots of the polynomial equation (12)

$$\begin{array}{rl}{\mathrm{\Delta }}_{1,1}=& {\displaystyle \frac{1}{a_4}(a_3{a}_4-a_2)>0,{\mathrm{\Delta }}_{1,2}={\displaystyle \frac{1}{a_4}(a_1{a}_4-a_0)>0,}}\\ {\mathrm{\Delta }}_{2,1}=& {\displaystyle \frac{1}{{\mathrm{\Delta }}_{1,1}}(a_2{\mathrm{\Delta }}_{1,1}-a_4{\mathrm{\Delta }}_{1,2})>0,{\mathrm{\Delta }}_{2,2}=a_0>0,{\mathrm{\Delta }}_{3,1}={\displaystyle \frac{1}{{\mathrm{\Delta }}_{2,1}}({\mathrm{\Delta }}_{1,2}{\mathrm{\Delta }}_{2,1}-a_0{\mathrm{\Delta }}_{1,1})>0.}}\end{array}$$

Consequently, the endemic equilibrium $\mathcal{E}*$ is locally asymptotically stable whenever ${\mathcal{R}}_{ef}>1.$

Uniform persistence

For the proposed model, the meningitis disease endemicity occurs if asymptomatic carriers, symptomatic infected and people undergoing medication fraction of the population persists above a certain positive level for sufficiently large time. This phenomenon can be well captured and analyzed through the notion of uniform persistence. And, system (1) is said to be uniformly persistent if there exists a constant $0<\xi <1$ such that any solution $(S,I_a,I_s,M,R)$ with initial conditions $(S^0,I_a^0,I_s^0,M^0,R^0)\in \mathrm{\Gamma }$ satisfies

$$\underset{t\to \mathrm{\infty }}{lim}infS(t)>\xi ,\underset{t\to \mathrm{\infty }}{lim}inf{I}_a(t)>\xi ,\underset{t\to \mathrm{\infty }}{\lim }inf{I}_s(t)>\xi ,\underset{t\to \mathrm{\infty }}{lim}infM(t)>\xi ,\underset{t\to \mathrm{\infty }}{lim}infR(t)>\xi$$

(13)

Consider the subsystem of the infected states of system (1) as follows:

$$\{\begin{array}{c}{\displaystyle \frac{d{I}_a}{dt}={\beta }_s{e}^{-m{I}_s}S{I}_s+{\beta }_{a}S{I}_a-(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )I_a,}\\ {\displaystyle \frac{d{I}_s}{dt}=(1-\gamma u_1)\alpha I_a-({\delta }_s+u_1+\mu )I_s,}\\ {\displaystyle \frac{dM}{dt}=\gamma u_1{I}_a+u_1{I}_s-({\delta }_h+{\gamma }_1{u}_2+\mu )M,}\end{array}$$

(14)and

$${\mathrm{\Gamma }}_1=\left\{(S,I_a,I_s,M,R):S>0,I_a\ge 0,I_s\ge 0,M\ge 0,R\ge 0,S+I_a+I_s+M+R<{\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\right\}\subset \mathrm{\Gamma },$$

where $\mathrm{\partial }{\mathrm{\Gamma }}_1$ the boundary of ${\mathrm{\Gamma }}_1.$ Let ${\mathrm{\Phi }}_t(x)=\mathrm{\Phi }(t,x(t))$ be the flow on ${\mathrm{\Gamma }}_1$ generated by the solution $x(t)$ of model (14) with the initial condition $x(0)=(I_a^0,I_s^0,M^0)\in {\mathrm{\Gamma }}_1.$

Theorem 3.5.

If ${\mathcal{R}}_{ef}>1$ , then the flow ${\mathrm{\Phi }}_t(x)$ on ${\mathrm{\Gamma }}_1$ is uniformly persistent for any Solutions $x(t)$ with the initial condition $I_a^0>0,I_s^0>0,M^0>0.$

Proof. We stated in Theorem 3.1 that the disease-free equilibrium is unstable when ${\mathcal{R}}_{ef}>1.$ Besides, one previously established that $\underset{t\to \mathrm{\infty }}{lim}supS(t)=\frac{\mathrm{\Lambda }}{\mu }$ , then there exists $ϵ>0$ such that

$$S(t)\le {\displaystyle \frac{\mathrm{\Lambda }}{\mu }+ϵ.}$$

(15)

Consequently, from model (14) we deduce the following auxiliary system:

$$\{\begin{array}{cc}& {\displaystyle \frac{d\widetilde{I_a}}{dt}={\beta }_s{e}^{-m\widetilde{I_s}}S\widetilde{I_s}+{\beta }_{a}S\widetilde{I_a}-(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )\widetilde{I_a},}\\ & {\displaystyle \frac{d\widetilde{I_s}}{dt}=(1-\gamma u_1)\alpha \widetilde{I_a}-({\delta }_s+u_1+\mu )\widetilde{I_s},}\\ & {\displaystyle \frac{d\tilde{M}}{dt}=\gamma u_1\widetilde{I_a}+u_1\widetilde{I_s}-({\delta }_h+{\gamma }_1{u}_2+\mu )\tilde{M},}\end{array}$$

(16)

Let

$$\mathrm{\Sigma }=\bigcup _{yϵY}\omega (y)$$

where

$$\omega (y)=\bigcap _{t\ge 0}{\overline{\mathcal{J}}}_t$$

with

$${\mathcal{J}}_t=\mathrm{\Phi }([0,+\mathrm{\infty }[)\times \{y\},t\in [0,+\mathrm{\infty }[$$

and

$$Y=\{x=(S,I_a,I_s,M,R):{\mathrm{\Phi }}_t(x)\in \mathrm{\partial }{\mathrm{\Gamma }}_1\}.$$

One can see that $\mathrm{\Sigma }$ is the maximal invariant set of the flow ${\mathrm{\Phi }}_t(x)$ on $\mathrm{\partial }{\mathrm{\Gamma }}_1$ . By analyzing the model system (1), we obtain that $\mathrm{\Sigma }$ consists of a unique equilibrium $dfe=\left(\frac{\mathrm{\Lambda }}{\mu },0,0,0,0\right)$ on the boundary $\mathrm{\partial }{\mathrm{\Gamma }}_1.$ Hence, the infection-free $dfe$ represents an acyclic covering of the set $\mathrm{\Sigma }.$

Now, let us analyze the behavior of any solution $x(t)$ of model (1) close to the disease-free point $dfe.$ For that, we classify the initial conditions in two cases:

If $I_a^0=I_s^0=M^0=0$ then $I_a(t)=I_s(t)=M(t)=0.$ Therefore, from System (16), we can see that any solution $x(t)$ of model (1) tends to $dfe.$

If $I_a^0>0,I_s^0>0$ and $M^0>0$ then by Proposition 3.1 $I_a(t)\ge 0,I_s\ge 0$ and $M(t)\ge 0.$ If the solution $x(t)$ is close to $dfe$ , according to system (1), there exists $\zeta >0$ such that

$$\{\begin{array}{c}{\displaystyle \frac{d{I}_a}{dt}>-(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu +\zeta )I_a,}\\ {\displaystyle \frac{d{I}_s}{dt}>[(1-\gamma u_1)\alpha -\zeta ]I_a-({\delta }_s+u_1+\mu +\zeta )I_s,}\\ {\displaystyle \frac{dM}{dt}>(\gamma u_1-\zeta )I_a+(u_1-\zeta )I_s-({\delta }_h+{\gamma }_1{u}_2+\mu +\zeta )M.}\end{array}$$

(17)

Therefore, system (17) reads as

$${\displaystyle \frac{d\mathcal{S}(t)}{dt}>\mathcal{B}\mathcal{S}(t)}$$

(18)where $\mathcal{S}=(I_a,I_s,M){}^T$ and

$$\mathcal{B}=(\begin{array}{ccc}-(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu +\zeta )& 0& 0\\ (1-\gamma u_1)\alpha -\zeta & -({\delta }_s+u_1+\mu +\zeta )& 0\\ \gamma u_1-\zeta & u_1-\zeta & -({\delta }_h+{\gamma }_1{u}_2+\mu +\zeta )\end{array}).$$

However, the largest eigenvalue of the matrix B is positive since ${\mathcal{R}}_{ef}>1$ .^4,15,28 Hence, the solutions of the quasi-monotonic linear system

$${\displaystyle \frac{d\tilde{\mathcal{S}}(t)}{dt}>\mathcal{B}\tilde{\mathcal{S}}(t)}$$

(19)with positive initial conditions are exponentially increasing as $t\to \mathrm{\infty }.$ By applying the comparison principle, we can see that $(I_a,I_s,M)$ diverges from $(0,0,0)$ when $t\to \mathrm{\infty }.$

Therefore, $\{dfe\}$ is an isolated invariant set of the flow ${\mathrm{\Phi }}_t(x).$ Thanks to the Theorem 4.3 in the work by Freedman et al., ³⁰ model (1) is uniformly persistent if ${\mathcal{R}}_{ef}>1$ . This completes the proof.

Global stability analysis of endemic equilibrium point

In this subsection, we investigate the global stability of the endemic equilibrium point.

Theorem 3.6. The unique endemic equilibrium

$\mathcal{E}*$ of model (1) is globally asymptotically stable when ${\mathcal{R}}_{ef}>1$ .

Proof.

Let us make use of the following Lyapunov candidate:

$$\mathcal{L}={\displaystyle \frac{1}{2}[(S-S*)+(I_a-I_a*)+(I_s-I_s*)+(M-M*)+(R-R*){]}^2}$$

(20)

It is clear that for all $(S,I_a,I_s,M,R)\in ({\mathbb{R}}_{+}^5)*,\mathcal{L}(S,I_a,I_s,M,R)>0$ and at the endemic equilibrium point, $\mathcal{L}(S*,I_a*,I_s*,M*,R*)=0.$

Now, by differentiating $\mathcal{L}$ with respect to time t alongside the solution path of the model (1), we obtain:

$$\begin{array}{rl}{\displaystyle \frac{d\mathcal{L}}{dt}=}& [(S-S*)+(I_a-I_a*)+(I_s-I_s*)+(M-M*)+(R-R*)]\times {\displaystyle \frac{d(S+I_a+I_s+M+R)}{dt}}\\ =& [(S+I_a+I_s+M+R)-(S*+I_a*+I_s*+M*+R*)]\times [\mathrm{\Lambda }-\mu (S+I_a+I_s+M+R)-{\delta }_s{I}_s-{\delta }_{h}M]\\ & \le [(S+I_a+I_s+M+R)-(S*+I_a*+I_s*+M*+R*)]\times [\mathrm{\Lambda }-\mu (S+I_a+I_s+M+R)]\\ & =-\mu [(S-S*)+(I_a-I_a*)+(I_s-I_s*)+(M-M*)+(R-R*)]\\ & \times \left[(S+I_a+I_s+M+R)-{\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\right]\le -\mu [(S-S*)+(I_a-I_a*)+(I_s-I_s*)+(M-M*)\\ & +(R-R*){]}^2,\left\{since-{\displaystyle \frac{\mathrm{\Lambda }}{\mu }\le -(S*+I_a*+I_s*+M*+R*)}\right\}\\ & {\displaystyle \frac{d\mathcal{L}}{dt}\le 0.}\end{array}$$

However, one can easily notice that $\frac{d\mathcal{L}}{dt}$ is semi-negative definite. Thus, $\mathcal{L}$ is a Lyapunov Function. Consequently, $\frac{d\mathcal{L}}{dt}\le 0$ with equality only at $\mathcal{E}*=(S*,I_a*,I_s*,M*,R*).$ It then follows by LaSalle's invariance principle that the largest invariant set $\mathrm{\Omega }$ , such that $\frac{d\mathcal{L}}{dt}=0$ is the singleton $\{\mathcal{E}*\}$ , this implies that the unique endemic equilibrium point is globally asymptotically stable.

Numerical results and sensitivity analysis

This section of the manuscript is devoted to the numerical simulation of model system (1) for both disease eradication (when ${\mathcal{R}}_0<1$ and ${\mathcal{R}}_{ef}<1$ ) and its persistence (when ${\mathcal{R}}_0>1$ and ${\mathcal{R}}_{ef}>1$ ) and local sensitivity analysis of the threshold parameter ${\mathcal{R}}_{ef}^a$ and ${\mathcal{R}}_{ef}^s$ that compose ${\mathcal{R}}_{ef}$ . In Table 2, we give the values of the parameters used to illustrate the cases of disease extinction and persistence. To begin, we will first examine the scenario in which no control measures are implemented, specifically when $u_1=u_2=0$ (see Figures 3 and 4). Next, we analyzed the case where control measures are applied at fixed levels, with $u_1=55\text{\%}$ and $u_2=45\text{\%}$ (see Figures 5 and 6).

Figure 3.

Simulations of model (1) in the case of disease extinction and without control measures $(u_1=u_2=0)$ ; when ${\mathcal{R}}_0\simeq 0.8567<1.$

Figure 4.

Simulations of model (1) in the case of disease persistence and without control measures $(u_1=u_2=0)$ ; when ${\mathcal{R}}_0\simeq 5.4759>1.$

Figure 5.

Simulations of model (1) in the case of disease extinction, with control measures $(u_1=0.55$ and $u_2=0.45)$ ; when ${\mathcal{R}}_{ef}\simeq 0.2708<1.$

Figure 6.

Simulations of model (1) in the case of disease persistence, with control measures $(u_1=0.55$ and $u_2=0.45)$ and ${\mathcal{R}}_{ef}\simeq 2.4298>1.$

Table 2.

Model parameter values.

Parameters	Values (extinction)	Values (persistence)	References
$\mathrm{\Lambda }$	360	450	27
${\beta }_s$	0.000065	0.00065	Assumed
${\beta }_a$	0.000085	0.00085	Assumed
$\mu$	0.1450	0.1450	4,27
${\delta }_s$	0.0195	0.195	4,23
${\delta }_h$	0.0159	0.8159	21,27
$\gamma$	0.85	0.9	Assumed
${\gamma }_1$	0.245	0.0245	Assumed
$\omega$	0.142	0.4142	21,25
$\alpha$	0.28	0.928	23
$\vartheta$	0.4	0.4	25
$m$	2	2	31

Model simulation results

Figures 3 and 5 display the evolution of the susceptible individuals, asymptomatic individuals, symptomatic infectious, hospitalized people and recovered ones when ${\mathcal{R}}_0\simeq 0.8567<1$ and ${\mathcal{R}}_{ef}\simeq 0.2708<1$ , respectively. One can observe that in both cases susceptible population tends to $\frac{\mathrm{\Lambda }}{\mu }$ while all the remaining state variables die out. This corroborates Theorem 3.2. However, when control measures are being implemented (Figure 5) the compartments $I_a$ , $I_s$ , M, and R populations will rapidly die out compared to the case where there is no control (Figure 3).

Figures 4 and 6 show the bacterial meningitis disease persistence within the population for ${\mathcal{R}}_0\simeq 5.4759>1$ and ${\mathcal{R}}_{ef}\simeq 2.4298>1$ , respectively. By considering the time-level up to 60 units (days), both susceptible individuals, asymptomatic individuals, symptomatic infectious and recovered persist within the population which corroborates Theorem 3.6. Moreover, Figures 5(d) and 6(d) represent two interesting situations: Figure 6(d) shows that in absence of rapid diagnosis and treatment, the hospitalized population will decrease to zero as time goes. However in the scenario (Figure 6(d)), where rapid diagnosis and treatment are both implemented up to 55% and 45%, respectively, the population of hospitalized patients will grow and reach some stability, which is well in keeping with Theorem 3.6.

Figure 7 displays the dynamics of asymptomatic people $I_a$ , symptomatic infectious $I_s$ , and isolated people M, for $m=0$ and $m=2.$ Clearly, the psychological impact of media coverage on people has an important impact in terms of slowing down the propagation of the disease transmission dynamics.

Figure 7.

Simulations of model (1) in the case of disease persistence and use of control measures $(u_1=0.55$ and $u_2=0.45)$ ; with and without media coverage influence on contacts of symptomatic infectious persons.

Sensitivity analysis of ${\mathcal{R}}_{ef}^a$ and ${\mathcal{R}}_{ef}^s$

The sensitivity analysis examines how variations in model inputs influence the model's outputs. Typically, two types of sensitivity analyses are conducted: local and global sensitivity analysis.

In this subsection, we focus on local sensitivity analysis, specifically through the lens of the effective reproduction number. For a parameter (p) that appears in both ${\mathcal{R}}_{ef}^a$ and ${\mathcal{R}}_{ef}^s$ , we compute the sensitivity index using Chitnis's normalized forward sensitivity index ⁴ as

$${\chi }_p^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^k}={\displaystyle \frac{\mathrm{\partial }{\mathcal{R}}_{ef}^k}{\mathrm{\partial }p}\times {\displaystyle \frac{p}{{\mathcal{R}}_{ef}^k}\simeq {\displaystyle \frac{\text{\%}\mathrm{\Delta }{\mathcal{R}}_{ef}^k}{\text{\%}\mathrm{\Delta }\mathrm{p}},k=a,s}}}$$

Figure 8 represents the local sensitivity indices from Tables 3 and 4 of the threshold parameters ${\mathcal{R}}_{ef}^a$ and ${\mathcal{R}}_{ef}^s$ that are the contribution of asymptomatic and symptomatic infected individuals, respectively. One can see that the parameters that appear within ${\mathcal{R}}_{ef}^a$ and ${\mathcal{R}}_{ef}^s$ differently impact the disease dynamics. Therefore, increasing $\mathrm{\Lambda },{\beta }_a$ or ${\beta }_s$ will speed the disease spreading while an increase in $\omega ,\alpha ,\gamma$ and $u_1$ values will slow down the disease progression.

Figure 8.

Local sensitivity diagram of ${\mathcal{R}}_{ef}^s$ (case (a)) and ${\mathcal{R}}_{ef}^a$ (case (b)).

Table 3.

Elasticity of ${\mathcal{R}}_{\mathit{e}\mathit{f}}^s$ with respect to various parameters.

Parameter	Symbol	Formula	Sensitivity index
$\mathrm{\Lambda }$	${\chi }_{\mathrm{\Lambda }}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^s}$	$-$	1
${\beta }_s$	${\chi }_{{\beta }_s}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^s}$	$-$	1
$\omega$	${\chi }_{\omega }^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^s}$	$-\frac{{\beta }_s\mathrm{\Lambda }(1-\gamma u_1)\alpha }{\mu ({\delta }_s+u_1+\mu ){\mathcal{Q}}^2}$	$-0.4580$
$\gamma$	${\chi }_{\gamma }^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^s}$	$-\frac{{\beta }_s\mathrm{\Lambda }(1-\gamma u_1)(1-\alpha )\alpha u_1}{\mu ({\delta }_s+u_1+\mu ){\mathcal{Q}}^2}$	$-0.0181$
$u_1$	${\chi }_{u_1}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^s}$	$-\frac{{\beta }_s\mathrm{\Lambda }(1-\gamma u_1)\alpha [\gamma (1-\alpha )({\delta }_s+u_1+\mu )+\mathcal{Q}}{\mu {({\delta }_s+u_1+\mu )}^2{\mathcal{Q}}^2}$	$-0.1611$

Table 4.

Elasticity of ${\mathcal{R}}_{\mathit{e}\mathit{f}}^a$ with respect to various parameters.

Parameter	Symbol	Formula	Sensitivity index
$\mathrm{\Lambda }$	${\chi }_{\mathrm{\Lambda }}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^a}$	$-$	1
${\beta }_a$	${\chi }_{{\beta }_a}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^a}$	$-$	1
$\omega$	${\chi }_{\omega }^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^a}$	$-\frac{{\beta }_a\mathrm{\Lambda }}{\mu {\mathcal{Q}}^2}$	$-1.1375$
$\alpha$	${\chi }_{\alpha }^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^a}$	$-\frac{{\beta }_a\mathrm{\Lambda }(1-\gamma u_1)}{\mu {\mathcal{Q}}^2}$	$-0.5744$
$\gamma$	${\chi }_{\gamma }^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^a}$	$-\frac{{\beta }_s\mathrm{\Lambda }(1-\alpha )u_1}{\mu {\mathcal{Q}}^2}$	$-0.0450$
$u_1$	${\chi }_{u_1}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^a}$	$-\frac{{\beta }_s\mathrm{\Lambda }(1-\alpha )\gamma }{\mu {\mathcal{Q}}^2}$	$-0.0737$

with

$$\mathcal{Q}=\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu .$$

By computing the sensitivity indices, we get for instance ${\chi }_{\omega }^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^s}=-0.4580$ and ${\chi }_{\omega }^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^a}=-1.1375$ that is $1\text{\%}$ increase in $\omega$ leads to a 0.4580% decrease in ${\mathcal{R}}_{ef}^s$ and a 1.1375% decrease in ${\mathcal{R}}_{ef}^a.$ Also, ${\chi }_{u_1}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^s}=-0.1611$ and ${\chi }_{u_1}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^a}=-0.0737$ mean that a $1\text{\%}$ increase in $u_1$ leads to a $0.1611\text{\%}$ decrease in ${\mathcal{R}}_{ef}^s$ and a $0.0737\text{\%}$ decrease in ${\mathcal{R}}_{ef}^a$ . Similarly ${\chi }_{{\beta }_s}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^s}=1$ , which means that a $1\text{\%}$ increase in ${\beta }_s$ only leads to a $1\text{\%}$ increase in ${\mathcal{R}}_{ef}^s$ , and ${\chi }_{{\beta }_a}^{{\mathcal{R}}_{\mathit{e}\mathit{f}}^a}=1$ means that a $1\text{\%}$ increase in ${\beta }_a$ only leads to a $1\text{\%}$ increase in ${\mathcal{R}}_{ef}^a.$ Similar interpretation can be done for the other parameters appearing in the expression of ${\mathcal{R}}_{ef}^s$ and ${\mathcal{R}}_{ef}^a.$

Figures 9 and 10 show the results of numerical simulations of model (1) for different values of $u_1$ and $\omega$ respectively, in situations where the disease persists. We observe that the higher the value of $u_1$ , the more asymptomatic carriers and symptomatic infectious populations decrease and the number of hospitalized individuals increases. Also, the higher the value of ω, the lower the number of asymptomatic carriers, symptomatic infectious individuals, and hospitalized individuals.

Figure 9.

Simulations of model (1) in the case of disease persistence, with $u_2=0.45$ and varying the value of $u_1$ .

Figure 10.

Simulations of model (1) in the case of disease persistence, with $u_1=0.55,u_2=0.45$ and varying the value of $\omega$ .

Optimal control analysis of the model when ${\mathcal{R}}_{ef}>1$

In this section, our purpose is to investigate the dynamical behavior of the proposed model (1) using Pontryagin's maximum principle ³² with a view to stemming the infectious population and keeping the requested cost of implementing early diagnosis $u_1(t)$ and treatment of hospitalized people $u_2(t)$ to the best bearable minimum. Otherwise, given that $y(t)\in {\mathrm{\Omega }}_y\subset {\mathbb{R}}^n$ is the state variable of model (1) and $u(t)\in {\mathrm{\Omega }}_u\subset {\mathbb{R}}^n$ are the control variables at any time t such that $0\le t\le T_f$ , then an optimal control problem consists of finding a piecewise continuous control $u(t)$ and its corresponding state $y(t).$

Furthermore, the Lagrangian of the system (1) is given by

$$\mathcal{L}(t,y,u)=A_1{I}_a+A_2{I}_s+A_{3}M+{\displaystyle \frac{1}{2}\sum _{i=1}^2{B}_i{u}_i^2}$$

(21)and the objective functional required to achieving this aim reads as

$$\mathcal{J}(u_1,u_2)=\underset{0}{\overset{T_f}{\int }}\left(A_1{I}_a+A_2{I}_s+A_{3}M+{\displaystyle \frac{1}{2}{B}_1{u}_1^2+{\displaystyle \frac{1}{2}{B}_2{u}_2^2}}\right)dt$$

(22)where $A_i,i=1,2,3$ are the balancing weight constant for the infectious individuals, $B_i,i=1,2$ are the weight constants for applying the control, and $T_f$ is the extent time frame for the control implementation. Further, the quantity $\underset{0}{\overset{T_f}{\int }}(A_1{I}_a+A_2{I}_s+A_{3}M)dt$ represents the cumulative cost due to the disease and the quantity $\underset{0}{\overset{T_f}{\int }}\frac{1}{2}(B_1{u}_1^2+B_2{u}_2^2)dt$ is the cumulative cost due to both early diagnosis and treatment of the infected people.

Therefore, using Pontryagin's Maximum Principle^32,33 this optimizes the functional $\mathcal{J}(u_1,u_2).$ Then, we seek optimal controls $u_1*$ and $u_2*$ that minimize the functional as follows:

$$\mathcal{J}(u_1*,u_2*)=\underset{(u_1,u_2)\in {\mathrm{\Omega }}_u}{min}\mathcal{J}(u_1,u_2)$$

(23)where Ω _u represents the Lebesgue measurable control set defined as

$${\mathrm{\Omega }}_u=\{(u_1,u_2):0\le u_i\le 1,t\in [0,T_f]\}$$

(24)

It is important to note that, the condition $u_i(t)=0$ implies that no control effort is being implemented while the condition $u_i(t)=1$ is the maximal control effort that is needed to successfully curtain the disease dynamics. These two cases could somewhat be too idealistic to implement.

Existence of the optimal control

Assuming the existence of a non-trivial vector $\lambda =({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_5)$ of adjoint or co-state variables satisfying the following conditions:

$${\displaystyle \frac{dy}{dt}={\displaystyle \frac{\mathrm{\partial }\mathbb{H}(t,y,u,\lambda )}{\mathrm{\partial }\lambda },{\displaystyle \frac{\mathrm{\partial }\mathbb{H}(t,y,u,\lambda )}{\mathrm{\partial }u}=0,{\displaystyle \frac{d\lambda }{dt}=-{\displaystyle \frac{\mathrm{\partial }\mathbb{H}(t,y,u,\lambda )}{\mathrm{\partial }y}}}}}}$$

(25)we achieve our goal by considering the Hamiltonian H for the control problem (1) expressed as

$$\begin{array}{rl}\mathbb{H}=& A_1{I}_a+A_2{I}_s+A_{3}M+{\displaystyle \frac{1}{2}{B}_1{u}_1^2+{\displaystyle \frac{1}{2}{B}_2{u}_2^2}}\\ & +{\lambda }_1[\mathrm{\Lambda }+\vartheta R-{\beta }_s{e}^{-m{I}_s}S{I}_s-{\beta }_{a}S{I}_a-\mu S]\\ & +{\lambda }_2[{\beta }_s{e}^{-m{I}_s}S{I}_s+{\beta }_{a}S{I}_a-(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )I_a]\\ & +{\lambda }_3[(1-\gamma u_1)\alpha I_a-({\delta }_s+u_1+\mu )I_s]\\ & +{\lambda }_4[\gamma u_1{I}_a+u_1{I}_s-({\delta }_h+{\gamma }_1{u}_2+\mu )M]\\ & +{\lambda }_5[{\gamma }_1{u}_{2}M+\omega I_a-(\vartheta +\mu )R].\end{array}$$

Now, by using equation (25) we derive the necessary conditions to the Hamiltonian H that our control functions and the corresponding state variables must satisfy.

Theorem 5.1. Consider the objective functional

$\mathcal{J}$ (22) defined on the control set ${\mathrm{\Omega }}_u$ subject to the optimal control system (1) with non-negative initial data at $t=0$ , there exists an optimal control $u*=(u_1*,u_2*)\in {\mathrm{\Omega }}_u$ satisfying equation (23).

Proof.

Given the control set ${\mathrm{\Omega }}_u=[0,1{]}^2,u=(u_1,u_2)\in {\mathrm{\Omega }}_u,x=(S*,I_a*,I_s*,M*,R*)$ and

$$f(t,x,u)=(\begin{array}{c}\mathrm{\Lambda }+\vartheta R-{\beta }_s{e}^{-m{I}_s}S{I}_s-{\beta }_{a}S{I}_a-\mu S\\ {\beta }_s{e}^{-m{I}_s}S{I}_s+{\beta }_{a}S{I}_a-(\omega +(1-\gamma u_1)\alpha +\gamma u_1+\mu )I_a\\ (1-\gamma u_1)\alpha I_a-({\delta }_s+u_1+\mu )I_s\\ \gamma u_1{I}_a+u_1{I}_s-({\delta }_h+{\gamma }_1{u}_2+\mu )M\\ {\gamma }_1{u}_{2}M+\omega I_a-(\vartheta +\mu )R\end{array})$$

be the right-hand side of our meningitis model (1).

Thanks to Fleming and Rishel Theorem, ³⁴ Theorem 5.1 hold if the following properties are preserved: $({\mathcal{A}}_1):$ The set of control ${\mathrm{\Omega }}_u$ is closed and convex, $({\mathcal{A}}_2):$ The model system (1) is bounded by a linear function in both the state and control variables, $({\mathcal{A}}_3):$ The Lagrangian of the objective functional (22) with respect to the control $(u_1*,u_2*)$ is convex and $({\mathcal{A}}_4):$ There exists constants $h_1,h_2>0$ and $h_3>1$ such that the Lagrangian is bounded below by $h_1({|u_i|}^2{)}^{h_3/2}-h_2,i=1,2.$

To establish the above properties $({\mathcal{A}}_1),({\mathcal{A}}_2),({\mathcal{A}}_3),({\mathcal{A}}_4)$ we proceed as follows:

$({\mathcal{A}}_1):$ Looking at the definition of the set ${\mathrm{\Omega }}_u=[0,1{]}^2$ of permissible controls, it is obvious to see that it is a closed set. Now, consider $p,q\in {\mathrm{\Omega }}_u$ two arbitrary points such that $p=(p_1,p_2)$ and $q=(q_1,q_2).$ Thanks to the definition of the convex set ex-rayed, ³⁴ given by

$$(\xi p_i+(1-\xi )q_i)\in [0,1],\mathrm{\forall }\xi \in [0,1],i=1,2.$$

Consequently, $(\xi p+(1-\xi )q)\in {\mathrm{\Omega }}_u$ implying that the control set ${\mathrm{\Omega }}_u$ is convex. Thereafter, property $({\mathcal{A}}_1)$ holds.

$({\mathcal{A}}_2):$ The proposed non-autonomous model (1) can be rewritten as a linear function of double control $u=(u_1,u_2)$ with respect to time and the coefficients of the dependent variables of the state system. This can only be achieved by using the explicit approach designed in Refs.^33,35 coupled with the boundedness of the solutions of the state system as stated by Proposition 3.1. Consequently, the expression on the right-hand side of the model (1) is bounded above by a sum of bounded state and control.

$({\mathcal{A}}_3):$ Here, we investigate the convexity of the Lagrangian of the objective functional (22) given by equation (21) with respect to the controls. It is given as

$$\mathcal{L}(t,y,u)=A_1{I}_a+A_2{I}_s+A_{3}M+{\displaystyle \frac{1}{2}\sum _{i=1}^2{B}_i{u}_i^2}$$

(26)where $y=(I_a,I_s,M)$ ,and $u=(u_1,u_2)$ . Therefore, for all $p=(p_1,p_2)\in {\mathrm{\Omega }}_u$ , $q=(q_1,q_2)\in {\mathrm{\Omega }}_u$ and $\xi \in [0,1]$ , it then follows that,

$$\mathcal{L}(t,y,\xi p+(1-\xi )q)=A_1{I}_a+A_2{I}_s+A_{3}M+{\displaystyle \frac{1}{2}\sum _{i=1}^2{B}_i(\xi p_i+(1-\xi )q_i{)}^2}$$

and

$$\xi \mathcal{L}(t,y,p)+(1-\xi )\mathcal{L}(t,y,q)=A_1{I}_a+A_2{I}_s+A_{3}M+{\displaystyle \frac{1}{2}\xi \sum _{i=1}^2{B}_i{p}_i^2+{\displaystyle \frac{1}{2}(1-\xi )\sum _{i=1}^2{B}_i{q}_i^2}}$$

Consequently, that leads to

$$\begin{array}{rl}& \mathcal{L}(t,y,\xi p+(1-\xi )q)-[\xi \mathcal{L}(t,y,p)+(1-\xi )\mathcal{L}(t,y,q)]\\ & ={\displaystyle \frac{1}{2}\sum _{i=1}^2{B}_i(\xi p_i+(1-\xi )q_i{)}^2-{\displaystyle \frac{1}{2}\xi \sum _{i=1}^2{B}_i{p}_i^2+{\displaystyle \frac{1}{2}(1-\xi )\sum _{i=1}^2{B}_i{q}_i^2}}}\\ & =-{\displaystyle \frac{1}{2}\xi (1-\xi )\sum _{i=1}^2{B}_i(p_i-q_i{)}^2}\\ & \le 0.\end{array}$$

Thus,

$$\mathcal{L}(t,y,\xi p+(1-\xi )q)\le \xi \mathcal{L}(t,y,p)+(1-\xi )\mathcal{L}(t,y,q)$$

(27)

And thereafter the Lagrangian is convex.

$({\mathcal{A}}_4):$ This property can easily be checked throughout the expression of the Lagrangian given by equation (21) by letting $h_1=min\{B_1/2,B_2/2\},h_2>0$ and $h_3=2.$

This completes the proof of the existence of an optimal control with respect to our optimal control problem.

Characterization of the optimal control

In this subsection, we derive the necessary conditions for the characterization needed to minimizing the objective functional (22) using Pontryagin's maximal principle. ³² Considering the Hamiltonian of the model (1), the following result yields

Theorem 5.2.

Given an optimal control measures $u_1*,u_2*$ and the steady states $S*,I_a*,I_s*,M*,R*$ of the corresponding state system (1), there exists adjoint variables $\lambda =({\lambda }_i{)}^T\in {\mathbb{R}}^5$ preserving the co-state system given by

$$\{\begin{array}{c}{\displaystyle \frac{d{\lambda }_1}{dt}=\mu {\lambda }_1+({\beta }_s{e}^{-m{I}_s}{I}_s+{\beta }_a{I}_a)({\lambda }_1-{\lambda }_2)}\\ {\displaystyle \frac{d{\lambda }_2}{dt}=-A_1+\mu {\lambda }_2+{\beta }_{a}S({\lambda }_1-{\lambda }_2)+(1-\gamma u_1)\alpha ({\lambda }_2-{\lambda }_3)+\gamma u_1({\lambda }_2-{\lambda }_4)+\omega ({\lambda }_2-{\lambda }_5)}\\ {\displaystyle \frac{d{\lambda }_3}{dt}=-A_2+({\delta }_s+\mu ){\lambda }_3+{\beta }_s(1-m{I}_s)e^{-m{I}_s}S({\lambda }_1-{\lambda }_2)+u_1({\lambda }_3-{\lambda }_4)}\\ {\displaystyle \frac{d{\lambda }_4}{dt}=-A_3+({\delta }_h+\mu ){\lambda }_4+{\gamma }_1{u}_2({\lambda }_4-{\lambda }_5)}\\ {\displaystyle \frac{d{\lambda }_5}{dt}=\mu {\lambda }_5+\vartheta ({\lambda }_5-{\lambda }_1)}\end{array}$$

(28)with the transversality conditions ${\lambda }_i(T_f)=0,i=1,2,\dots ,5$ and the optimal control characterized as follows:

$$u_1*=max\left\{min\left\{{\displaystyle \frac{({\lambda }_3-{\lambda }_2)\alpha \gamma I_a*+({\lambda }_2-{\lambda }_4)\gamma I_a*+({\lambda }_3-{\lambda }_4)I_s*}{B_1},u_1^{max}}\right\},0\right\}$$

(29)

$$u_2*=max\left\{min\left\{{\displaystyle \frac{({\lambda }_4-{\lambda }_5){\gamma }_{1}M*}{B_2},u_2^{max}}\right\},0\right\}$$

(30)

Proof. Indeed, differentiating the Hamiltonian of the system (1) with respect to the variables of the state system leads to

$${\displaystyle \frac{d{\lambda }_1}{dt}=-{\displaystyle \frac{\mathrm{\partial }\mathbb{H}}{\mathrm{\partial }S},{\displaystyle \frac{d{\lambda }_2}{dt}=-{\displaystyle \frac{\mathrm{\partial }\mathbb{H}}{\mathrm{\partial }{I}_a},{\displaystyle \frac{d{\lambda }_3}{dt}=-{\displaystyle \frac{\mathrm{\partial }\mathbb{H}}{\mathrm{\partial }{I}_s},{\displaystyle \frac{d{\lambda }_4}{dt}=-{\displaystyle \frac{\mathrm{\partial }\mathbb{H}}{\mathrm{\partial }M},{\displaystyle \frac{d{\lambda }_5}{dt}=-{\displaystyle \frac{\mathrm{\partial }\mathbb{H}}{\mathrm{\partial }R}.}}}}}}}}}}$$

Moreover, to characterize the optimal control one solves through the following optimality conditions for $u_i*,i=1,2$

$${\displaystyle \frac{\mathrm{\partial }\mathbb{H}}{\mathrm{\partial }{u}_1}{|}_{u_1*}=0and{\displaystyle \frac{\mathrm{\partial }\mathbb{H}}{\mathrm{\partial }{u}_2}{|}_{u_2*}=0.}}$$

In addition, each control $u_i*,i=1,2$ satisfies these standard control arguments involving the limits of the controls

$$u_i*=\{\begin{array}{cc}0& if{\phi }_i*\le 0,\\ {\phi }_i*& if0<{\phi }_i*\le u_i^{max},\\ u_i^{max}& if{\phi }_i*\ge u_i^{max},\end{array}$$

for $i=1,2$ and with

$${\phi }_1*={\displaystyle \frac{({\lambda }_3-{\lambda }_2)\alpha \gamma I_a*+({\lambda }_2-{\lambda }_4)\gamma I_a*+({\lambda }_3-{\lambda }_4)I_s*}{B_1}and{\phi }_2*={\displaystyle \frac{({\lambda }_4-{\lambda }_5){\gamma }_{1}M*}{B_2}.}}$$

Optimal control simulation

Here, we perform the numerical simulation of the optimal control model by employing the forward-backward iterative method. Besides, we consider a situation where the disease is persistent that is ${\mathcal{R}}_0>1.$ And, the weight values of the objective function (22) are given by $A_1=1.2,A_2=1.2,A_3=1.5,B_1=1,B_2=1$ and the initial conditions are given as follows: $S(0)=S^0=1215,I_a(0)=I_a^0=170,I_s(0)=I_s^0=15,\mathrm{M}(0)=M^0=150,\mathrm{R}(0)=R^0=147.$ We also used the disease persistence parameter values given in Table 2. As stated before, two control measures are considered namely the control $u_1(t)$ for early diagnosis and $u_2(t)$ patients management in hospitals. Considering the importance of both of these two control strategies and their complementary effects, we display the impact of implementing them both together. We used the software MATLAB and the fourth-order Runge-Kutta iterative scheme to solve the optimal system (1)–(22) and the equations satisfied by the adjoint variables, in order to obtain the optimal control solution. First, the forward difference approximation is applied to solve the optimal system (1), then the backward difference approximation is used to solve the system (28), of equations satisfied by the adjoint variables, using the solutions of the current iterations of the state equations and considering the transversality conditions. Then, at each iteration, we update the control functions using the optimality conditions given in (29) and (30), then the iterations continue until convergence.^{23,24,33–36} For the readers convenience, we remind that we marked the dynamical behavior of the model without control by red lines and the dynamics with control are marked by blue lines. Figure 11 clearly shows that implementing some control strategies can slow the disease spreading. For instance, by the early diagnosis, asymptomatic and symptomatic infectious can be isolated for medication and prevent them sharing the disease (see Figures 11(b) to (d)). Moreover, the adjoint variables and the optimal control profiles are given in Figures 12 and 13, respectively.

Figure 11.

Simulations of model (1) in the case of application of optimal control strategies ( $u_1$ and $u_2$ ) to a situation of disease persistence $({\mathcal{R}}_0=5.4759)$ .

Figure 12.

Adjoints variables.

Figure 13.

Control profiles.

Discussion and conclusion

In conclusion, this comprehensive mathematical analysis highlights the critical role of various factors, including media coverage, early diagnosis, isolation, and treatment, in preventing the bacterial meningitis epidemic. Through the application of sophisticated methodologies such as the next-generation method to computing the effective reproduction number, the Castillo- Chavez theorem for the global stability of the infection-free point when ${\mathcal{R}}_{ef}<1$ , the Lyapunov functional technique for the global stability if ${\mathcal{R}}_{ef}>1$ , and the maximum Pontryagin principle to investigate the optimal problem, we have developed a robust and multifaceted framework that enhances our understanding of the complex dynamics involved in meningitis disease transmission and control.^4,6,28,33

Our results reveal that timely and effective communication strategies are paramount to raise public awareness and encourage proactive health-seeking behaviors. The role of the media in disseminating accurate information cannot be overemphasized, as they help inform the public about the symptoms of bacterial meningitis, the importance of early diagnosis, and the steps to take in case of suspected infection. In parallel, rapid implementation of medical interventions, including isolation of affected individuals and appropriate treatment protocols, is essential to curb the spread of this serious infection.^11,12,14

In addition, the local sensitivity analysis shows that the parameters whose variations induce more change in the value of ${\mathcal{R}}_{ef}$ are: the recruitment rate $\mathrm{\Lambda }$ , the contact rate between susceptible individuals and asymptomatic carriers (resp. symptomatic individuals) ${\beta }_a$ (resp. ${\beta }_s$ ), the natural recovery rate of asymptomatic carriers $\omega$ , the rate $\alpha$ at which asymptomatic carriers become symptomatic infectious individuals, followed by the early diagnosis rate $u_1$ . For this constructed model, the local sensitivity study then suggests that measures to combat bacterial meningitis should aim to reduce the values of the parameters $\mathrm{\Lambda },{\beta }_a,{\beta }_s$ and increase those of the parameters $\omega ,u_1,\alpha$ , and $\gamma$ . Figures 9 and 10 are perfect illustrations of these conclusions. We also illustrated (Figure 7) that information relayed by the media can contribute to reducing bacterial meningitis infections if followed by the population. The study of optimal control and its numerical simulation reveal the possibility of curbing the spread of bacterial meningitis through early diagnosis and patient care, while minimizing the cost of implementation. However, to achieve these results, a synergy of actions would be required so that well-informed populations can be diagnosed as early as possible and infectious individuals can be isolated and treated for the benefit of all. This study provides a better understanding of the factors that influence the epidemic dynamics of meningitis. By modeling the interactions between the different components of disease spread, we can identify critical thresholds and optimal intervention strategies that public health decision-makers can employ. This analysis not only contributes to our theoretical understanding of bacterial meningitis epidemics, but also serves as a practical guide for developing targeted strategies to preserve community health.^4,6

In summary, our study emphasizes the interdependence of media influence, healthcare accessibility, and rapid intervention in the fight against bacterial meningitis. By leveraging mathematical frameworks, we can better equip public health officials with the knowledge and tools needed to prevent future outbreaks, ultimately leading to better health outcomes and increased community resilience to infectious diseases.^15,20,21

Looking ahead, this model can serve as a foundation for developing a bacterial meningitis model tailored to any specific country, situated within the meningitis belt incorporating factors like seasonality, spatial dynamics, vaccination, and age-structured. Using real data, it will enable to perform accurate simulations to inform recommendations that support refining the national strategy and achieving the target meningitis eradication by 2030.