Mathematical Modeling and Analysis of Ebola Virus Disease Dynamics: Implications for Intervention Strategies and Healthcare Resource Optimization

Abstract

This study implements a minded approach to studying Ebola virus disease (EVD) by dividing the infected population into aware and unaware groups and including a hospitalized compartment. This offers a more detailed understanding of illness distribution, potential analyses, and the influence of public knowledge. The findings might improve healthcare budget apportionment, public health policy, and contest Ebola and related infections. In this study, we fully observe the new model $SEIHR$ that we have constructed. We start by outlining the essential concepts of the model and confirming its mathematical reliability. Next, we calculate the fundamental reproductive number $\left(R_0\right)$, which is critical for appreciating how the infection spreads and how effective treatments might be. We also study stability analysis, which looks at when the disease may decline or become chronic. Furthermore, we exhibit the occurrence of bifurcation in the EVD Epidemic Model and perform a sensitivity analysis of $\left(R_0\right)$. The main findings of this study show that for $R_0<1$, the disease-free equilibrium, is globally stable, meaning the disease will die out, whereas for $R_0>1$, the endemic equilibrium is stable, meaning the disease persists. Additionally, the sensitivity analysis reveals that the most influential parameters in controlling $R_0$ are the transmission rate and the recovery rate, which could guide effective intervention strategies. Finally, we use numerical simulations so that out outcomes are more significant.

1. Introduction

Ebola virus disease (EVD) is caused by a virus that most likely originated in flying foxes, often known as fruit bats. These bats can transfer and spread the virus to humans [1]. The first cases of Ebola were noted in 1976, following severe fever epidemics in Central Africa. Interacting with animals such as porcupines, monkeys, apes, and fruit bats, which are natural hosts of the Ebola virus, increases the risk of transmission. These animals may transfer the virus without symptoms, serving as reservoirs in their natural environments. When people come into contact with them, particularly in areas where Ebola is prevalent, they risk contracting the virus. Contact with biological fluids such as blood, saliva, and urine, as well as activities such as hunting, can spread the infection. Understanding these transmission routes is essential for preventing outbreaks and protecting public health. Educating people about the risks of handling these animals and encouraging safe practices in doing so is important for containing the spread of Ebola and reducing its impact on both human and animal communities. Healthcare workers and aid staff can also contract the virus while treating patients in healthcare facilities [2,3]. Health authorities have launched workshops and educational programs for doctors in Pakistan [4], acknowledging the epidemic threat of Ebola.

Various fields have used mathematical modeling to observe different aspects of diseases and their dynamics [5,6,7,8,9,10,11,12,13,14,15]. These models can encompass a wide range of medical resources or focus on specific elements such as hospital beds and laboratory equipment. For example, Zhou and Fan [5] examined human behavior in the presence of scarce medical resources using an $SIR$ model, which takes into account susceptibility, infection, and recovery. Their results revealed that individuals tend to adopt more protective measures when faced with limited treatment and prevention options for highly deadly diseases.

The Ebola virus disease has been widely studied using mathematical models that consider limited resources. In a study by Nouvellet et al. [6], an $SEIR$ model was developed, incorporating delays between different stages of the disease, to highlight the importance of rapid diagnostic tests for Ebola. Similarly, Leward et al. [7] examined pharmacological interventions to control EVD and arrived at similar conclusions. In 2020, Rafiq et al. [14] analyzed a model that classified the population into four groups: susceptible individuals ($S\left(t\right)$), exposed individuals ($E\left(t\right)$), infected individuals ($I\left(t\right)$), and individuals who have recovered ($R\left(t\right)$). Their focus was on evaluating the stability of the different equilibrium points of the model [14]. Previous studies have explored various strategies for controlling Ebola, emphasizing the importance of proactive measures in managing outbreaks [10,11,16].

This study introduces key innovations compared to prior work. Specifically, we divide the infected population into two distinct groups: “aware” and “unaware” individuals, based on their knowledge of the disease. Furthermore, we include a hospitalized compartment to account for the role of healthcare facilities in controlling the outbreak. Unlike previous models, which primarily focused on general population compartments, our model provides a more granular analysis of public health interventions. This approach allows us to assess the impact of public awareness campaigns and hospital capacity on disease dynamics, making our contribution unique in highlighting the importance of behavioral and structural interventions in managing EVD outbreaks.

The rest of the paper is organized as follows: Section 2 outlines the model formulation. Then, Section 3 discusses how the model is well-posed, with subsections focusing on specific aspects like Invariant Regions and the Positivity of Solutions. Section 4 discusses Steady States, with subsections on the Disease-Free Equilibrium Point (DFE) and the Endemic Equilibrium Point (EE). The Basic Reproductive Number ($R_0$) is computed in Section 5. Section 6 and Section 7 analyze the stability of the model, both locally and globally, considering equilibrium points. Section 8 discusses the existence of Bifurcation within the model. Sensitivity Analysis of $R_0$ is performed in Section 9, followed by Numerical Simulations in Section 10. Finally, Section 11 concludes the paper.

2. Model Formulation

In order to understand how Ebola spreads between humans and animals, we have modified $SEIHR$ model for Ebola transmission. This model categorizes the total population into six groups over time. $S\left(t\right)$ represents susceptible individuals who are at risk of contracting Ebola, $E\left(t\right)$ denotes the exposed class who have been in contact with the Ebola virus but are not yet showing symptoms, $I_u\left(t\right)$ signifies infected individuals who are unaware of their infection at the current time, $I_a\left(t\right)$ indicates infected individuals who are aware of their infection at the current time, $H\left(t\right)$ represents hospitalized individuals who require medical care, and $R\left(t\right)$ denotes recovered individuals who are no longer infectious. The flowchart and parameter descriptions are given in Figure 1 and

Table 1. Parameter descriptions.

Parameter	Description
$\Lambda$	New recruitment rate
$\mu$	Natural mortality rate
${\delta }_1$	Death rate due to disease for exposed individuals
${\delta }_2$	Mortality rate due to disease for unaware individuals
${\delta }_3$	Mortality rate due to disease for aware individuals
${\delta }_4$	Mortality rate due to disease for hospitalized individuals
${\alpha }_1$	Rate of infection from susceptible to exposed humans
${\alpha }_2$	Transmission rate from exposed to unaware persons
${\alpha }_3$	Rate of transmission from exposed to aware peoples
${\alpha }_4$	Transmission rate from individuals in the unaware state to those in the hospitalized state
${\alpha }_5$	Transmission rate from individuals in the aware state to those in the hospitalized state
${\alpha }_6$	Transmission rate from individuals in the unaware state to those in the recovered state
${\alpha }_7$	Transmission rate from individuals in the aware state to those in the recovered state
${\alpha }_8$	Transmission rate from individuals in the hospitalized state to those in the recovered state
${\eta }_1$	Infection rate from wild animals to individuals in the exposed state
${\eta }_2$	Infection rate from wildlife to individuals in the unaware state
${\eta }_3$	Infection rate from wildlife to individuals in the aware state
${\sigma }_1$	Infection rate from domestic animals to individuals in the exposed state
${\sigma }_2$	Infection rate from domestic animals to individuals in the unaware state
${\sigma }_3$	Infection rate from domestic animals to individuals in the aware state

, respectively. Through this model, we track the movement of individuals between these categories to understand the dynamics of Ebola transmission among humans and animals. In this study, the following model is considered:

$$\left\{\begin{array}{c}{\displaystyle \frac{dS}{dt}}=\Lambda -\mu S-({\alpha }_1+{\eta }_1+{\sigma }_1)SE-({\sigma }_2+{\eta }_2)S{I}_u-({\sigma }_3+{\eta }_3)S{I}_a,\hfill \\ {\displaystyle \frac{dE}{dt}}=({\alpha }_1+{\eta }_1+{\sigma }_1)SE-{\alpha }_2{I}_{u}E-{\alpha }_3{I}_{a}E-(\mu +{\delta }_1)E,\hfill \\ {\displaystyle \frac{d{I}_u}{dt}}=({\sigma }_2+{\eta }_2)S{I}_u+{\alpha }_2{I}_{u}E-{\alpha }_{4}H{I}_u-(\mu +{\delta }_2+{\alpha }_6)I_u,\hfill \\ {\displaystyle \frac{d{I}_a}{dt}}=({\sigma }_3+{\eta }_3)S{I}_a+{\alpha }_3{I}_{a}E-{\alpha }_{5}H{I}_a-(\mu +{\delta }_3+{\alpha }_7)I_a,\hfill \\ {\displaystyle \frac{dH}{dt}}={\alpha }_{4}H{I}_a+{\alpha }_{5}H{I}_u-(\mu +{\delta }_4+{\alpha }_8)H,\hfill \\ {\displaystyle \frac{dR}{dt}}={\alpha }_6{I}_u+{\alpha }_7{I}_a+{\alpha }_{8}H-\mu R.\hfill \end{array}\right.$$

(1)

with initial conditions: $S\left(0\right)>0$, and $(E\left(0\right),I_u\left(0\right),I_a\left(0\right),H\left(0\right),R\left(0\right))\ge 0$. The mathematical model is illustrated in the flowchart in Figure 1, and the description of the parameters is provided in Table 1.

3. Well-Posedness of the Proposed Model

3.1. Invariant Region

Theorem 1.

Within the region $\tau =\{(S,E,I_u,I_a,H,R)\in {\mathbb{R}}_{+}^6:S>0,(E,I_u,I_a,H,R)\ge 0,N\le {\displaystyle \frac{\Lambda }{\mu }}\}$, the system described by Equation (1) exhibits the existence and uniqueness of a solution given specific initial conditions. This region τ remains both attracting and positively invariant concerning the system’s dynamics. The proof of its existence and uniqueness can be established using Picard’s existence theorem.

Proof. 

The functions on the right-hand side of (1) are Lipschitz continuous, and they play a crucial role in Picard’s existence theorem, which determines when solutions to the system are possible. This means that the system described by (1) has solutions that are both positive and bounded within a certain region $\tau$. Now, let us talk about the total population, denoted by N and is given by

$$N=S+E+I_u+I_a+H+R.$$

(2)

By differentiating Equation (2) with respect to time and then putting the values from (1), we obtain

$${\displaystyle \frac{dN}{dt}}+\mu N\le \Lambda .$$

(3)

Solving (3) through the application of an integrating factor and incorporating the initial condition $N\left(0\right)=N_0$, we obtain

$$N\le N_0{e}^{-\mu t}+{\displaystyle \frac{\Lambda }{\mu }}(1-e^{-\mu t}).$$

(4)

In (4), as t approaches infinity, N tends to ${\displaystyle \frac{\Lambda }{\mu }}$. Consequently, we observe that in system (1), the entire set of feasible solutions converges to the enclosed region $\tau$. Thus, the closed region $\tau$ remains positively invariant at all times t, affirming the biological significance and mathematical coherence of model (1) within $\tau$.  □

3.2. Positivity of the Solutions

Our model aims to ensure biological realism. We aim to show that every path starting from the non-negative region ${\mathbb{R}}_{+}^6$ will ultimately converge to and stay within the feasible area $\tau$. To prove this, we will establish that the set $\tau$ is positively invariant and is the system’s global attractor [17].

Theorem 2.

Consider thagt the initial conditions of (1) are nonnegative, then the solutions for the different groups $S,E,I_u,I_a,H$ and R of system (1) remain nonnegative $\forall t>0$.

Proof. 

Take into consideration, $t^{\prime }=sup\{t>0;S>0,$ and $(E,I_u,I_a,H,R)\ge 0\}$.

In that case, $t^{\prime }>0,$ because $S\left(0\right)>0$, also take

$$\begin{array}{c}\hfill {\displaystyle \frac{dS}{dt}}=\Lambda -S[\mu +({\sigma }_1+{\eta }_1+{\alpha }_1)E+({\sigma }_2+{\eta }_2)I_u+({\sigma }_3+{\eta }_3)I_a].\end{array}$$

(5)

Let $g\left(t\right)=({\sigma }_1+{\eta }_1+{\alpha }_1)E+({\sigma }_2+{\eta }_2)I_u+({\sigma }_3+{\eta }_3)I_a$, then (5) becomes

$$\begin{array}{c}\hfill {\displaystyle \frac{dS}{dt}}+(\mu +g\left(t\right))S=\Lambda .\end{array}$$

(6)

After solving (6) and using $S\left(0\right)=S_0$, yields

$$\begin{array}{c}\hfill S\left(t^{\prime }\right)=S_0{e}^{(-\mu t^{\prime }-{\int }_0^{t^{\prime }}g\left(t\right)dt)}+\Lambda e^{(-\mu t^{\prime }-{\int }_0^{t^{\prime }}f\left(t\right)dt)}{\int }_0^{t^{\prime }}{e}^{(\mu y+{\int }_0^{y}g\left(x\right)dx)}dy.\end{array}$$

(7)

As $S_0>0$, $\Lambda >0$ also $e^{(-\mu t^{\prime }-{\int }_0^{t^{\prime }}g\left(t\right)dt)}$, $e^{(-\mu t^{\prime }-{\int }_0^{t^{\prime }}f\left(t\right)dt)}$, and ${\int }_0^{t^{\prime }}{e}^{(\mu y+{\int }_0^{y}g\left(x\right)dx)}dy$ are positive. Therefore, $S\left(t\right)>0$.

From model (1), Equation (2) can be written as

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dE}{dt}}& \ge & -E({\alpha }_2{I}_u+{\alpha }_3{I}_a)-(\mu +{\delta }_1)E.\hfill \end{array}$$

(8)

Let $g_1\left(t\right)={\alpha }_2{I}_a+{\alpha }_3{I}_u$. By solving (8) and using $E\left(0\right)=E_0$, we obtain

$$E\ge E_0{e}^{(\mu +{\mu }_1)t}{e}^{\int -\left(g_1\left(t\right)\right)dt}.$$

(9)

In (9) $E_0\ge 0$, $e^{\int -\left(g_1\left(t\right)\right)dt}>0$ and $e^{(\mu +{\mu }_1)t}>0.$ and Therefore $E\ge 0$.

A similar procedure can be extended to the rest of the equations in model (1). This analysis shows that all the state variables we are tracking stay non-negative for any $t\ge 0$. This is important because it means our model’s solutions always remain positive, which is crucial for its accuracy.  □

4. Equilibrium Points

In disease modeling, steady states refer to situations where the number of people in each group remains constant over time. This happens when the rates at which people move between groups are balanced. It means the number of infected individuals remains stable over time for a specific value of $R_0$, no matter how many people were initially infected. Here, we discuss two types of steady states (disease-free equilibrium and endemic equilibrium).

4.1. Disease-Free Equilibrium Point (DFE)

To determine the DFE $P_0$ of System (1), each equation on the right-hand side of system (1) is set to zero. After simplifying and arranging the expressions, we obtain $P_0=(S_0,E_0,I_{u_0},I_{a_0},H_0,R_0)=({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0,0)$.

4.2. Endemic Equilibrium Point (EE)

Similarly, to determine the EE point $P_1=(S_1,E_1,I_{u_1},I_{a_1},H_1,R_1)$ of System (1), each equation on the right-hand side of system (1) is set as zero. Simplifying and arranging formulas for the endemic point $P_1$ yields the following result:

$$S_1={\displaystyle \frac{\Lambda ({\alpha }_4{\alpha }_3-{\alpha }_5{\alpha }_2)}{{\alpha }_4(\mu {\alpha }_3+b_1{b}_6-b_3{b}_4)+{\alpha }_5(b_4{b}_2-b_1{b}_5-\mu {\alpha }_2)+b_7({\alpha }_3{b}_2-{\alpha }_2{b}_3)}},$$

$$\begin{array}{l}{E}_1={\displaystyle \frac{\begin{array}{cc}\hfill & {\alpha }_4^2({\alpha }_3(\mu b_6-\Lambda b_3)+b_1{b}_6^2-b_3{b}_4{b}_6)+{\alpha }_4({\alpha }_5({\alpha }_2(\Lambda b_3-\mu b_6)+{\alpha }_3(\Lambda b_2-\mu b_5)+b_5(b_3{b}_4-2b_1{b}_6)+b_2{b}_4{b}_6)\hfill \\ \hfill & +b_6{b}_7({\alpha }_3{b}_2-{\alpha }_2{b}_3))+{\alpha }_5({\alpha }_5({\alpha }_2(\mu b_5-\Lambda b_2)+b_1{b}_5^2-b_2{b}_4{b}_5)-b_5{b}_7({\alpha }_3{b}_2-{\alpha }_2{b}_3))\hfill \end{array}}{({\alpha }_4(\mu {\alpha }_3+b_1{b}_6-b_3{b}_4)+{\alpha }_5(b_4{b}_2-b_1{b}_5-\mu {\alpha }_2)+b_7({\alpha }_3{b}_2-{\alpha }_2{b}_3))({\alpha }_3{\alpha }_4-{\alpha }_2{\alpha }_5)}},\hfill \\ \\ I_{u_1}={\displaystyle \frac{\begin{array}{cc}\hfill & {\alpha }_5^2(b_4{b}_2^2-b_4(b_5{b}_1+\mu {\alpha }_2)+{\alpha }_2\Lambda b_1)+{\alpha }_5({\alpha }_3({\alpha }_4(\mu b_4-\Lambda b_1)-b_7(\mu {\alpha }_2+b_5{b}_1)-2b_4{b}_2)\hfill \\ \hfill & +b_4\left({\alpha }_4(b_6{b}_1-b_4{b}_3-{\alpha }_2{b}_7{b}_3)\right)+{\alpha }_3{b}_7({\alpha }_3(\mu {\alpha }_4+b_7{b}_2)+{\alpha }_4(b_6{b}_1-b_4{b}_3)-{\alpha }_2{b}_7{b}_3)\hfill \end{array}}{\begin{array}{cc}\hfill & ({\alpha }_5(b_4{b}_2-b_5{b}_1-\mu {\alpha }_2)+{\alpha }_3(\mu {\alpha }_4+b_7{b}_2)+{\alpha }_4(b_6{b}_1-b_4{b}_3)-{\alpha }_2{b}_7{b}_3)({\alpha }_3{\alpha }_4-{\alpha }_2{\alpha }_5)\hfill \end{array}}},\hfill \\ \\ I_{a_1}={\displaystyle \frac{\begin{array}{cc}\hfill & {\alpha }_4^2({\alpha }_3(\Lambda b_1-\mu b_4)-b_4{b}_1{b}_6+b_3{b}_4^2)+{\alpha }_4({\alpha }_2({\alpha }_5(\mu b_4-\Lambda b_1)-b_7(\mu {\alpha }_3+b_6{b}_1-2b_3{b}_4))\hfill \\ \hfill & +b_4\left({\alpha }_5(b_5{b}_1-b_2{b}_4)\right)-{\alpha }_3{b}_2{b}_7)+{\alpha }_2{b}_7({\alpha }_2(\mu {\alpha }_5+b_7{b}_3)+{\alpha }_5(b_5{b}_1-b_2{b}_4)-{\alpha }_3{b}_2{b}_7)\hfill \end{array}}{\begin{array}{cc}\hfill & ({\alpha }_4(\mu b_3+b_6{b}_1-b_4{b}_3)-{\alpha }_2(\mu {\alpha }_5+b_7{b}_3)+{\alpha }_5(b_4{b}_2-b_5{b}_1)+{\alpha }_3{b}_2{b}_7)({\alpha }_3{\alpha }_4-{\alpha }_2{\alpha }_5)\hfill \end{array}}},\hfill \\ \\ H_1={\displaystyle \frac{\begin{array}{cc}\hfill & {\alpha }_2^2({\alpha }_5(\Lambda b_3-\mu b_6)-b_3{b}_6{b}_7)+{\alpha }_2({\alpha }_3({\alpha }_4(\mu b_6-\Lambda b_3)+{\alpha }_5(\mu b_5-\Lambda b_2)+b_7(b_2{b}_6+b_3{b}_5))+b_6({\alpha }_4(b_1{b}_6-b_3{b}_4)\hfill \\ \hfill & -{\alpha }_5(b_1{b}_5-b_2{b}_4)\left)\right)-{\alpha }_3({\alpha }_3({\alpha }_4(\mu b_5-\Lambda b_2)+b_2{b}_5{b}_7)+b_5({\alpha }_4(b_1{b}_6-b_3{b}_4)-{\alpha }_5(b_1{b}_5-b_2{b}_4)))\hfill \end{array}}{\begin{array}{cc}\hfill & ({\alpha }_2(-\mu {\alpha }_5-b_3{b}_7)+{\alpha }_3(\mu {\alpha }_4+b_2{b}_7)+{\alpha }_4(b_1{b}_6-b_3{b}_4)-{\alpha }_5(b_1{b}_5-b_2{b}_4))({\alpha }_3{\alpha }_4-{\alpha }_2{\alpha }_5)\hfill \end{array}}},\hfill \\ \\ R_1={\displaystyle \frac{\begin{array}{cc}\hfill & {\alpha }_2^2({\alpha }_5(\mu {\alpha }_7{b}_7+{\alpha }_8(\Lambda b_3-\mu b_6))+b_3{b}_7({\alpha }_7{b}_7-{\alpha }_8{b}_6))+{\alpha }_4({\alpha }_3({\alpha }_4(-\mu {\alpha }_7{b}_7-{\alpha }_8(\Lambda b_3-\mu b_6))+{\alpha }_5(\mu {\alpha }_6{b}_7\hfill \\ \hfill & -{\alpha }_8(\Lambda b_2-\mu b_5))-b_7(b_7({\alpha }_6{b}_3+{\alpha }_7{b}_2)-{\alpha }_8(b_2{b}_6+b_3{b}_5)))+{\alpha }_5(-{\alpha }_7(\Lambda b_1-\mu b_4))+b_4(2{\alpha }_7{b}_3{b}_7-{\alpha }_8{b}_3{b}_6)\hfill \\ \hfill & -b_1{b}_6({\alpha }_7{b}_7-{\alpha }_8{b}_6))+{\alpha }_5({\alpha }_6(\Lambda b_1-\mu b_4){\alpha }_5+b_4(b_7(-{\alpha }_6{b}_3-{\alpha }_7{b}_2)+{\alpha }_8{\alpha }_2{b}_6)+b_1{b}_5({\alpha }_7{b}_7-{\alpha }_8{b}_6)){\alpha }_2\hfill \\ \hfill & +{\alpha }_3^2({\alpha }_4(\mu {\alpha }_6{b}_7+{\alpha }_8(\Lambda b_2-\mu b_5))+b_2{b}_7({\alpha }_6{b}_7-{\alpha }_8{b}_5))+{\alpha }_3({\alpha }_7(\Lambda b_1-\mu b_4){\alpha }_4^2+{\alpha }_4(-{\alpha }_6(\Lambda b_1-\mu b_4){\alpha }_5\hfill \\ \hfill & +b_4(b_7(-{\alpha }_6{b}_3-{\alpha }_7{b}_2)+{\alpha }_8{b}_3{b}_5)+b_1{b}_6({\alpha }_6{b}_7-{\alpha }_8{b}_5))-{\alpha }_5(b_4(-2{\alpha }_6{b}_2{b}_7+{\alpha }_8{b}_2{b}_5)+b_1{b}_5({\alpha }_6{b}_7-{\alpha }_8{b}_5)))\hfill \\ \hfill & -({\alpha }_4{\alpha }_7-{\alpha }_5{\alpha }_6){\alpha }_4({\alpha }_4(b_1{b}_6-b_3{b}_4)-{\alpha }_5(b_1{b}_5-b_2{b}_4))\hfill \end{array}}{\begin{array}{cc}\hfill & ({\alpha }_2(-\mu {\alpha }_5-b_3{b}_7)+{\alpha }_3(\mu {\alpha }_4+b_2{b}_7)+{\alpha }_4(b_1{b}_6-b_3{b}_4)-{\alpha }_5(b_1{b}_5-b_2{b}_4))\mu ({\alpha }_3{\alpha }_4-{\alpha }_2{\alpha }_5)\hfill \end{array}}},\hfill \end{array}$$

where $b_1=({\eta }_1+{\alpha }_1+{\sigma }_1)$, $b_2=({\eta }_2+{\sigma }_2)$, $b_3=({\eta }_3+{\sigma }_3)$, $b_4=(\mu +{\delta }_1)$, $b_5=(\mu +{\delta }_2+{\alpha }_6)$, $b_6=(\mu +{\delta }_3+{\alpha }_7)$, and $b_7=(\mu +{\delta }_4+{\alpha }_8)$.

5. Basic Reproduction Number $\left({\mathit{R}}_{\mathbf{0}}\right)$

To find the basic reproduction number $R_0$ for model (1), the next-generation matrix method is employed. From model (1), attention is directed towards the infected classes, following procedures outlined in various papers such as [18]. Matrices F and V are as follows.

$$F=\left[\begin{array}{ccc}& ({\eta }_1+{\alpha }_1+{\sigma }_1)SE-{\alpha }_2{I}_{u}E-{\alpha }_3{I}_{a}E& \\ & ({\eta }_2+{\sigma }_2)S{I}_u+{\alpha }_2{I}_{u}E-{\alpha }_{3}H{I}_u& \\ & ({\eta }_3+{\sigma }_3)S{I}_a+{\alpha }_3{I}_{a}E-{\alpha }_{5}H{I}_a& \\ & {\alpha }_{4}H{I}_u+{\alpha }_{5}H{I}_a\end{array}\right],V=\left[\begin{array}{ccc}& (\mu +{\delta }_1)E& \\ & (\mu +{\delta }_2+{\alpha }_6)I_u& \\ & (\mu +{\delta }_3+{\alpha }_7)I_a& \\ & (\mu +{\delta }_4+{\alpha }_8)H\end{array}\right].$$

The Jacobian matrices for F and V at the disease-free equilibrium point are presented as $F^{*}$ and $V^{*}$ are given below;

$$F^{*}=\left[\begin{array}{cccc}{\displaystyle \frac{\Lambda }{\mu }}({\eta }_1+{\sigma }_1+{\alpha }_1)& 0& 0& 0\\ 0& {\displaystyle \frac{\Lambda }{\mu }}({\eta }_2+{\sigma }_2)& 0& 0\\ 0& 0& {\displaystyle \frac{\Lambda }{\mu }}({\eta }_3+{\sigma }_3)& 0\\ 0& 0& 0& 0\end{array}\right],V^{*}=\left[\begin{array}{cccc}(\mu +{\delta }_1)& 0& 0& 0\\ 0& (\mu +{\delta }_2+{\alpha }_6)& 0& 0\\ 0& 0& (\mu +{\delta }_3+{\alpha }_7)& 0\\ 0& 0& 0& (\mu +{\delta }_4+{\alpha }_8)\end{array}\right].$$

Thus, the next-generation matrix is given by

$$F^{*}·V^{{*}^{-1}}=\left[\begin{array}{cccc}{\displaystyle \frac{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}{\mu (\mu +{\delta }_1)}}& 0& 0& 0\\ 0& {\displaystyle \frac{\Lambda ({\eta }_2+{\sigma }_2)}{\mu (\mu +{\delta }_2+{\alpha }_6)}}& 0& 0\\ 0& 0& {\displaystyle \frac{\Lambda ({\eta }_3+{\sigma }_3)}{\mu (\mu +{\delta }_3+{\alpha }_7)}}& 0\\ 0& 0& 0& 0\end{array}\right].$$

This matrix is diagonal, and thus its eigenvalues are as follows:

$$\begin{array}{c}\hfill {\lambda }_1={\displaystyle \frac{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}{\mu (\mu +{\delta }_1)}},{\lambda }_2={\displaystyle \frac{\Lambda ({\eta }_2+{\sigma }_2)}{\mu (\mu +{\delta }_2+{\alpha }_6)}},{\lambda }_3={\displaystyle \frac{\Lambda ({\eta }_3+{\sigma }_3)}{\mu (\mu +{\delta }_3+{\alpha }_7)}},{\lambda }_4=0.\end{array}$$

The basic reproduction number, designated as $R_0$, and derived as the dominant eigenvalue of $F^{*}·V^{{*}^{-1}}$. On completing the necessary computation, it is found that $R_0=max\{{\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4\}={\lambda }_1$, where ${\lambda }_i,i=1,2,3,4$ are eigenvalues derived from the calculation of the resulting matrix. Therefore,

$$R_0={\displaystyle \frac{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}{\mu (\mu +{\delta }_1)}}.$$

(10)

Graphical representations are currently illustrating how parameters influence the basic reproductive number $R_0$.

Figure 2a,b demonstrate that $R_0$ is directly proportional to $\Lambda$, ${\eta }_1$, ${\alpha }_1$, and ${\sigma }_1$. In Figure 2c, it is depicted that $R_0$ varies inversely with $\mu$ and ${\delta }_1$.

6. Local Stability Analysis of the Disease-Free Equilibrium (DFE) Point

This section entails conducting a local stability analysis of the DFE point $P_0$ for the model (1). We follow procedures outlined by various researchers [19,20] at $P_0=({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0,0)$, with some particular initial population values of each class.

Theorem 3.

The model (1) is locally asymptotically stable at $P_0$, within τ, when $R_0<1$, and unstable if $R_0>1$.

Proof. 

The Jacobian for system (1) at $P_0$ is below:

$$J\left(P_0\right)=\left(\begin{array}{cccccc}-\mu & -{\displaystyle \frac{\Lambda }{\mu }}{b}_1& -{\displaystyle \frac{\Lambda }{\mu }}{b}_2& -{\displaystyle \frac{\Lambda }{\mu }}{b}_3& 0& 0\\ 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_1-b_4& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_2-b_5& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_3-b_6& 0& 0\\ 0& 0& 0& 0& -b_7& 0\\ 0& 0& {\alpha }_6& {\alpha }_7& {\alpha }_8& -\mu \end{array}\right).$$

where $b_1=({\eta }_1+{\alpha }_1+{\sigma }_1)$, $b_2=({\eta }_2+{\sigma }_2)$, $b_3=({\eta }_3+{\sigma }_3)$, $b_4=(\mu +{\delta }_1)$, $b_5=(\mu +{\delta }_2+{\alpha }_6)$, $b_6=(\mu +{\delta }_3+{\alpha }_7)$, and $b_7=(\mu +{\delta }_4+{\alpha }_8)$.

From the above matrix, the eigenvalues are:

$$\begin{array}{c}{\lambda }_1=-\mu ,{\lambda }_2=-\mu ,{\lambda }_3=-(\mu +{\delta }_4+{\alpha }_8),{\lambda }_4={\displaystyle \frac{\Lambda }{\mu }}({\eta }_3+{\sigma }_3)-(\mu +{\delta }_3+{\alpha }_7),{\lambda }_5={\displaystyle \frac{\Lambda }{\mu }}({\eta }_2+{\sigma }_2)-(\mu +{\delta }_2+{\alpha }_6),\hfill \\ {\lambda }_6=(\mu +{\delta }_1)(R_0-1).\hfill \end{array}$$

Here, ${\lambda }_1=-\mu ,{\lambda }_2=-\mu ,$ because $\mu >0$, ${\lambda }_3=-(\mu +{\delta }_4+{\alpha }_8)<0$ as $(\mu +{\delta }_4+{\alpha }_8)>0$, ${\lambda }_4={\displaystyle \frac{\Lambda }{\mu }}({\eta }_3+{\sigma }_3)-(\mu +{\delta }_3+{\alpha }_7)<0$ if and only if ${\displaystyle \frac{\Lambda }{\mu }}({\eta }_3+{\sigma }_3)<(\mu +{\delta }_3+{\alpha }_7)$, ${\lambda }_5={\displaystyle \frac{\Lambda }{\mu }}({\eta }_2+{\sigma }_2)-(\mu +{\delta }_2+{\alpha }_6)<0\iff {\displaystyle \frac{\Lambda }{\mu }}({\eta }_2+{\sigma }_2)<(\mu +{\delta }_2+{\alpha }_6)$, and finally, if $R_0<1$, then all eigenvalues are less than 0, implying that $P_0=({\displaystyle \frac{\Lambda }{\mu }},0,0,0,0,0)$ is LAS, otherwise unstable.  □

7. Global Stability Analysis

This portion entails conducting global stability of system (1). For this, we use the Lyapunov function method to establish global asymptotic stability (GAS) at both equilibrium points, following the approach outlined in [21,22].

7.1. Global Stability at DFE Point

Here, the main goal is to construct the Lyapunov function, a specialized mathematical formula that ensures system stability and facilitates reaching the desired state even when conditions change.

Theorem 4.

The disease-free equilibrium point $P_0$ of the $EVD$ system (1) is globally asymptotically stable in τ when $R_0<1.$

Proof. 

Demonstrating the global stability of model (1) at $P_0$, we consider $V:\tau \to \mathbb{R}$ [20] as a Lyapunov function:

$$V=S-S_0(lnS)+E+I_u+I_a+H.$$

(11)

By differentiating (11) with respect to t and substituting the values from (1) into (11), we obtain:

$$\begin{array}{ccc}\hfill {\displaystyle \frac{dV}{dt}}& =& -{\displaystyle \frac{\mu }{S}}{(S-S_0)}^2-(\mu +{\delta }_4+{\alpha }_8)H+\left[(\mu +{\delta }_1)(R_0-1)\right]E\hfill \\ & +& [{\displaystyle \frac{\Lambda }{\mu }}({\eta }_2+{\sigma }_2)-(\mu +{\delta }_3+{\alpha }_7)]I_u+[{\displaystyle \frac{\Lambda }{\mu }}({\eta }_3+{\sigma }_3)-(\mu +{\delta }_2+{\alpha }_6)]I_a.\hfill \end{array}$$

(12)

As ${\displaystyle \frac{\Lambda }{\mu }}({\eta }_2+{\sigma }_2)<(\mu +{\delta }_3+{\alpha }_7)$ and ${\displaystyle \frac{\Lambda }{\mu }}({\eta }_2+{\sigma }_2)<(\mu +{\delta }_2+{\alpha }_6)$ from Theorem 3, ${\displaystyle \frac{dV}{dt}}\le 0$ for $R_0<1$. Definitely, ${\displaystyle \frac{dV}{dt}}=0$ only at $(S,E,I_u,I_a,H)=(S_0,0,0,0,0)$, and ${\displaystyle \frac{dV}{dt}}<0$ otherwise. This shows ${\displaystyle \frac{dV}{dt}}$ is negative semi-definite around $P_0=(S_0,0,0,0,0,0)$, approving V as a Lyapunov function in $\tau$. The point $P_0$ is the largest invariant subset where ${\displaystyle \frac{dV}{dt}}=0$.

Anywhere a solution trajectory that starts within the region $\tau$ with specific initial conditions, it gradually converges to $P_0$ as $t\to +\infty$. Consequently, the prevalence of Ebola virus disease (EVD) steadily decreases from the host population. Therefore, according to La Salle’s invariance principles [23], we can determine that $P_0$ exhibits global asymptotic stability within the domain $\tau$.  □

7.2. Global Stability Analysis of Endemic Equilibrium Point

Our objective here is to come up with a Lyapunov function that facilitates achieving global asymptotic stability.

Theorem 5.

If $R_0>1$, the equilibrium point $P_1$ of model (1) shows Global Asymptotic Stability in region τ for values $S_1$, $E_1$, $I_{u_1}$, $I_{a_1}$, and $H_1$. Conversely, if $R_0<1$, the model becomes unstable.

Proof. 

The Lyapunov function $V_1:\tau \to \mathbb{R}$ [20] is considered for showing the global stability of model (1) at the EE point $P_1$, defined as:

$$\begin{array}{cc}\hfill V_1& =K_1(S-S_{1}lnS)+K_2(E-E_{1}lnE)+K_3(I_u-I_{u_1}ln{I}_u)\hfill \\ \hfill & +K_4(I_a-I_{a_1}ln{I}_a)+K_5(H-H_{1}lnH),\hfill \end{array}$$

(13)

where $K_1$, $K_2$, $K_3$, $K_4$, and $K_5$ represent positive constants that will be determined later. When investigating the solution of model (1), we compute the rate of change of $V_1$ over time, giving the following expression:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{V}_1}{dt}}& =K_1(S-S_1)[-({\eta }_1+{\alpha }_1+{\sigma }_1)E-({\eta }_2+{\sigma }_2)I_u-({\eta }_3+{\sigma }_3)I_a-\mu ]\hfill \\ \hfill & +K_2(E-E_1)[({\eta }_1+{\alpha }_1+{\sigma }_1)S-{\alpha }_2{I}_u-{\alpha }_3{I}_a-(\mu +{\delta }_1)]\hfill \\ \hfill & +K_3(I_u-I_{u_1})[({\eta }_2+{\sigma }_2)S+{\alpha }_{2}E-{\alpha }_{4}H-(\mu +{\delta }_2+{\alpha }_6)]\hfill \\ \hfill & +K_4(I_a-I_{a_1})[({\eta }_3+{\sigma }_3)S+{\alpha }_{3}E-{\alpha }_{5}H-(\mu +{\delta }_3+{\alpha }_7)]\hfill \\ \hfill & +K_5(H-H_1)[{\alpha }_4{I}_u+{\alpha }_5{I}_a-(\mu +{\delta }_4+{\alpha }_8)].\hfill \end{array}$$

(14)

As $P_1=(S_1,E_1,I_{u_1},I_{a_1},H_1,R_1)$ represents an EE point, by substituting the LHS of all equations of (1) as follows:

$${\displaystyle \frac{d{S}_1}{dt}}={\displaystyle \frac{d{E}_1}{dt}}={\displaystyle \frac{d{I}_{u_1}}{dt}}={\displaystyle \frac{d{I}_{a_1}}{dt}}={\displaystyle \frac{dH}{dt}}=0,$$

which implies that:

$$\left\{\begin{array}{c}\mu ={\displaystyle \frac{\Lambda }{S_1}}-({\eta }_1+{\alpha }_1+{\sigma }_1)E_1-({\eta }_2+{\sigma }_2)I_{u_1}-({\eta }_3+{\sigma }_3)I_{a_1},\hfill \\ \mu +{\delta }_1=({\eta }_1+{\alpha }_1+{\sigma }_1)S_1-{\alpha }_2{I}_{u_1}-{\alpha }_3{I}_{a_1},\hfill \\ \mu +{\delta }_2+{\alpha }_6=({\eta }_2+{\sigma }_2)S_1+{\alpha }_2{E}_1-{\alpha }_4{H}_1,\hfill \\ \mu +{\delta }_3+{\alpha }_7=({\eta }_3+{\sigma }_3)S_1+{\alpha }_3{E}_1-{\alpha }_5{H}_1,\hfill \\ \mu +{\delta }_4+{\alpha }_8={\alpha }_4{I}_{u_1}+{\alpha }_5{I}_{a_1}.\hfill \end{array}\right.$$

(15)

Substituting (15) in (14) and taking $K_1=K_2=K_3=K_4=K_5=1$, we get:

$${\displaystyle \frac{d{V}_1}{dt}}=K_1(S-S_1)({\displaystyle \frac{\Lambda }{S}}-{\displaystyle \frac{\Lambda }{S_1}}).$$

This implies:

$${\displaystyle \frac{d{V}_1}{dt}}=-\Lambda {\displaystyle \frac{{(S-S_1)}^2}{S{S}_1}}\le 0.$$

Hence, $V_1$ represents a Lyapunov function. The condition ${\displaystyle \frac{d{V}_1}{dt}}=0$ remains valid if $(S,E,I_a,I_u,H,R)=(S_1,E_1,I_{u_1},I_{a_1},H_1,R_1)=P_1$, thus indicating that the point $P_1$ accounts for the most extensive invariant subset in set:

$${\epsilon }_1=\{(S,E,I_u,I_a,H,R)\in \tau :{\displaystyle \frac{d{V}_1}{dt}}=0\}.$$

This means that anywhere a solution trajectory that begins within the feasible region, $\tau$ gradually converges to $P_1$ as $t\to +\infty$. Consequently, the prevalence of EVD steadily increases from the host population. Hence, according to La Salle’s invariance principle [23], we can determine that $P_1$ exhibits global asymptotic stability within the domain $\tau$.  □

8. Bifurcation Existence in EVD Epidemic Model

Bifurcation in mathematical models involves small parameter changes causing notable changes in system behavior. It identifies critical points where system dynamics fundamentally change, aiding our understanding of system evolution under varying conditions.

For bifurcation analysis, we will apply central manifold theory, as detailed in [24]. To investigate bifurcation in our model, we set $R_0=1$ to determine $b_1^{*}$, where $b_1^{*}={\eta }_1+{\alpha }_1+{\sigma }_1$.

$$\begin{array}{ccc}\hfill R_0& =& {\displaystyle \frac{\Lambda b_1^{*}}{\mu (\mu +{\delta }_1)}}=1.\hfill \\ \hfill \Rightarrow b_1^{*}& =& {\displaystyle \frac{\mu (\mu +{\delta }_1)}{\Lambda }}={\displaystyle \frac{\mu b_4}{\Lambda }}.\hfill \end{array}$$

where $b_4=\mu +{\delta }_1$, the Jacobian matrix for model (1) at $P_0$ in DFE proceeds as below:

$$J^{\mid 0\mid }=\left(\begin{array}{cccccc}-\mu & -{\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}& -{\displaystyle \frac{\Lambda }{\mu }}{b}_2& -{\displaystyle \frac{\Lambda }{\mu }}{b}_3& 0& 0\\ 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}-b_4& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_2-b_5& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_3-b_6& 0& 0\\ 0& 0& 0& 0& -b_7& 0\\ 0& 0& {\alpha }_6& {\alpha }_7& {\alpha }_8& -\mu \end{array}\right).$$

Let us define the following parameters:

$$\begin{array}{cc}\hfill b_1& =({\eta }_1+{\alpha }_1+{\sigma }_1),b_2=({\eta }_2+{\sigma }_2),b_3=({\eta }_3+{\sigma }_3),b_4=(\mu +{\delta }_1),b_5=(\mu +{\delta }_2+{\alpha }_6),\hfill \\ \hfill b_6& =(\mu +{\delta }_3+{\alpha }_7),b_7=(\mu +{\delta }_4+{\alpha }_8).\hfill \end{array}$$

Upon calculation, the Jacobian matrix $J^{\mid 0\mid }$ eigenvalues are given as follows:

$$\begin{array}{cc}\hfill {\lambda }_{1,2}& =-\mu ,\hfill \\ \hfill {\lambda }_3& =-(\mu +{\delta }_4+{\alpha }_8),\hfill \\ \hfill {\lambda }_4& ={\displaystyle \frac{\Lambda }{\mu }}({\eta }_3+{\sigma }_3)-(\mu +{\delta }_3+{\alpha }_7),\hfill \\ \hfill {\lambda }_5& ={\displaystyle \frac{\Lambda }{\mu }}({\eta }_2+{\sigma }_2)-(\mu +{\delta }_2+{\alpha }_6),\hfill \\ \hfill {\lambda }_6& =(\mu +{\delta }_1)(R_0-1).\hfill \end{array}$$

For $R_0=1$, ${\lambda }_6=0$. Thus, we have one zero eigenvalue, indicating that a bifurcation exists.

Identifying the Nature of Bifurcation

This subsection analyzes the nature of the bifurcation presented in our model. To do so, the process involves identifying both the left and right eigenvectors, and then, based on these eigenvectors, ascertaining the type of bifurcation present, as outlined in [24].

$$J^{\left|0\right|}\left(b_1^{*}\right).U=0.$$

$$\left(\begin{array}{cccccc}-\mu & -{\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}& -{\displaystyle \frac{\Lambda }{\mu }}{b}_2& -{\displaystyle \frac{\Lambda }{\mu }}{b}_3& 0& 0\\ 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}-b_4& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_2-b_5& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_3-b_6& 0& 0\\ 0& 0& 0& 0& -b_7& 0\\ 0& 0& {\alpha }_6& {\alpha }_7& {\alpha }_8& -\mu \end{array}\right)\left(\begin{array}{c}{u}_1\\ u_2\\ u_3\\ u_4\\ u_5\\ u_6\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right),$$

which gives:

$$\left(\begin{array}{c}-\mu u_1-{\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}{u}_2-{\displaystyle \frac{\Lambda }{\mu }}{b}_2{u}_3-{\displaystyle \frac{\Lambda }{\mu }}{b}_3{u}_4\\ \left({\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}-b_4\right)u_2\\ \left({\displaystyle \frac{\Lambda }{\mu }}{b}_2-b_5\right)u_3\\ \left({\displaystyle \frac{\Lambda }{\mu }}{b}_3-b_6\right)u_4\\ -b_7{u}_5\\ {\alpha }_6{u}_3+{\alpha }_7{u}_4+{\alpha }_8{u}_5-\mu u_6\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right).$$

Therefore:

$$\begin{array}{c}-\mu u_1-{\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}{u}_2-{\displaystyle \frac{\Lambda }{\mu }}{b}_2{u}_3-{\displaystyle \frac{\Lambda }{\mu }}{b}_3{u}_4=0,\end{array}$$

(16)

$$\begin{array}{c}\left({\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}-b_4\right)u_2=0,\end{array}$$

(17)

$$\begin{array}{c}\left({\displaystyle \frac{\Lambda }{\mu }}{b}_2-b_5\right)u_3=0,\end{array}$$

(18)

$$\begin{array}{c}\left({\displaystyle \frac{\Lambda }{\mu }}{b}_3-b_6\right)u_4=0,\end{array}$$

(19)

$$\begin{array}{c}-b_7{u}_5=0,\end{array}$$

(20)

$$\begin{array}{c}{\alpha }_6{u}_3+{\alpha }_7{u}_4+{\alpha }_8{u}_5-\mu u_6=0.\end{array}$$

(21)

After solving Equations (16)–(21), we have $u_1=0$, $u_2=0$, $u_3=0$, $u_4=0$, $u_5=0$, and $u_6=0$. Therefore, the right eigenvector is zero.

For the zero eigenvalues, the left eigenvector is as follows:

$$C.J^{\left|0\right|}\left(b_1^{*}\right)=0.$$

$$\left(\begin{array}{cccccc}{c}_1& c_2& c_3& c_4& c_5& c_6\end{array}\right)\left(\begin{array}{cccccc}-\mu & -{\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}& -{\displaystyle \frac{\Lambda }{\mu }}{b}_2& -{\displaystyle \frac{\Lambda }{\mu }}{b}_3& 0& 0\\ 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}-b_4& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_2-b_5& 0& 0& 0\\ 0& 0& 0& {\displaystyle \frac{\Lambda }{\mu }}{b}_3-b_6& 0& 0\\ 0& 0& 0& 0& -b_7& 0\\ 0& 0& {\alpha }_6& {\alpha }_7& {\alpha }_8& -\mu \end{array}\right)=Q,$$

where

$$Q=\left(\begin{array}{cccccc}0& 0& 0& 0& 0& 0\end{array}\right).$$

After multiplying these matrices and equating them to zero, we obtain:

$$\begin{array}{c}-\mu c_1=0,\end{array}$$

(22)

$$\begin{array}{c}-{\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}{c}_1+\left({\displaystyle \frac{\Lambda }{\mu }}{b}_1^{*}-b_4\right)c_2=0,\end{array}$$

(23)

$$\begin{array}{c}-{\displaystyle \frac{\Lambda }{\mu }}{b}_2{c}_1+\left({\displaystyle \frac{\Lambda }{\mu }}{b}_2-b_5\right)c_3+{\alpha }_6{c}_6=0,\end{array}$$

(24)

$$\begin{array}{c}-{\displaystyle \frac{\Lambda }{\mu }}{b}_3{c}_1+\left({\displaystyle \frac{\Lambda }{\mu }}{b}_3-b_6\right)c_4+{\alpha }_7{c}_6=0,\end{array}$$

(25)

$$\begin{array}{c}-b_7{c}_5+{\alpha }_8{c}_6=0,\end{array}$$

(26)

$$\begin{array}{c}-\mu v_6=0.\end{array}$$

(27)

When we solve Equations (22)–(27), we have $c_1=0$, $c_2=0$, $c_3=0$, $c_4=0$, $c_5=0$, and $c_6=0$. Thus, the left eigenvector also equals zero. As both the right and left eigenvectors are zero, this represents a trivial case.

9. Sensitivity Analysis of ${\mathit{R}}_{\mathbf{0}}$

In this study, conducting a sensitivity analysis on model parameters regarding $R_0$ is crucial for identifying key factors influencing disease transmission and control. We evaluate parameter significance relative to $R_0$ using the formula by Chitnis et al. [25]. The definition of the sensitivity index on $R_0$ is

$$\begin{array}{c}{{\rm Y}}_{\kappa }^{R_0}={\displaystyle \frac{\partial R_0}{\partial \kappa }}\times {\displaystyle \frac{\kappa }{R_0}},\kappa \in \{\Lambda ,\mu ,{\delta }_1,{\eta }_1,{\alpha }_1,{\sigma }_1\}.\hfill \end{array}$$

$$\begin{array}{c}{{\rm Y}}_{\Lambda }^{R_0}={\displaystyle \frac{\partial }{\partial \Lambda }}\left({\displaystyle \frac{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}{\mu (\mu +{\delta }_1)}}\right)\times {\displaystyle \frac{\Lambda \mu (\mu +{\delta }_1)}{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}}=1,\hfill \\ {{\rm Y}}_{\mu }^{R_0}={\displaystyle \frac{\partial }{\partial \mu }}\left({\displaystyle \frac{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}{\mu (\mu +{\delta }_1)}}\right)\times {\displaystyle \frac{{\mu }^2(\mu +{\delta }_1)}{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}}=-1-{\displaystyle \frac{\mu }{\mu +{\delta }_1}}=-1.7854\hfill \\ {{\rm Y}}_{{\delta }_1}^{R_0}={\displaystyle \frac{\partial }{\partial {\delta }_1}}\left({\displaystyle \frac{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}{\mu (\mu +{\delta }_1)}}\right)\times {\displaystyle \frac{{\delta }_1\mu (\mu +{\delta }_1)}{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}}=-{\displaystyle \frac{{\delta }_1}{\mu +{\delta }_1}}=-0.2145,\hfill \\ {{\rm Y}}_{{\eta }_1}^{R_0}={\displaystyle \frac{\partial }{\partial {\eta }_1}}\left({\displaystyle \frac{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}{\mu (\mu +{\delta }_1)}}\right)\times {\displaystyle \frac{{\eta }_1\mu (\mu +{\delta }_1)}{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}}={\displaystyle \frac{{\eta }_1}{{\eta }_1+{\alpha }_1+{\sigma }_1}}=0.1521,\hfill \\ {{\rm Y}}_{{\alpha }_1}^{R_0}={\displaystyle \frac{\partial }{\partial {\alpha }_1}}\left({\displaystyle \frac{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}{\mu (\mu +{\delta }_1)}}\right)\times {\displaystyle \frac{{\alpha }_1\mu (\mu +{\delta }_1)}{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}}={\displaystyle \frac{{\alpha }_1}{{\eta }_1+{\alpha }_1+{\sigma }_1}}=0.3547,\hfill \\ {{\rm Y}}_{{\sigma }_1}^{R_0}={\displaystyle \frac{\partial }{\partial {\sigma }_1}}\left({\displaystyle \frac{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}{\mu (\mu +{\delta }_1)}}\right)\times {\displaystyle \frac{{\sigma }_1\mu (\mu +{\delta }_1)}{\Lambda ({\eta }_1+{\alpha }_1+{\sigma }_1)}}={\displaystyle \frac{{\sigma }_1}{{\eta }_1+{\alpha }_1+{\sigma }_1}}=0.4931.\hfill \end{array}$$

The sensitivity analysis of $R_0$ concerning the model parameters is illustrated in Figure 3. Sensitivity analysis of $R_0$ concerning all parameters encompasses both numerical computation and graphical representation. In these parameters, $\Lambda$, ${\eta }_1$, ${\alpha }_1$, and ${\sigma }_1$ exert a positive influence on $R_0$, whereas $\mu$ and ${\delta }_1$ exhibit a negative impact. Therefore, one can conclude that the most sensitive parameters are $\Lambda$ and $\mu$.

10. Numerical Simulations

This section focuses on the numerical simulations of the proposed model. To better understand how the model system (1) behaves dynamically in different scenarios, we conducted numerical simulations using the RK4 method and MATLAB R2023b. These simulations used the parameter values listed in

Table 2. Assigned values to the parameters of system (1) for simulation purposes.

Parameter	Value	Source	Parameter	Value	Source
$\Lambda$	$0.63210$	[14]	$\mu$	$0.43740$	Assumed
${\delta }_1$	$0.20060$	[14]	${\delta }_2$	$0.00210$	[26]
${\delta }_3$	$0.00130$	Assumed	${\delta }_4$	$0.00120$	Assumed
${\alpha }_1$	$0.28770$	[14]	${\alpha }_2$	$0.10000$	[26]
${\alpha }_3$	$0.00190$	[26]	${\alpha }_4$	$0.00100$	[26]
${\alpha }_5$	$0.00050$	[26]	${\alpha }_6$	$0.00210$	[26]
${\alpha }_7$	$0.00140$	Assumed	${\alpha }_8$	$0.00010$	[26]
${\eta }_1$	$0.12340$	[14]	${\eta }_2$	$0.24310$	[14]
${\eta }_3$	$0.08100$	Assumed	${\sigma }_1$	$0.40000$	[14]
${\sigma }_2$	$0.30000$	[14]	${\sigma }_3$	$0.10000$	Assumed

. In addition to time-series plots, we focused on revealing key dynamical behaviors such as bifurcation curves, limit cycles, and other relevant dynamical structures to provide a comprehensive understanding of the model’s dynamics. Appropriate parameter ranges were selected to visualize these phenomena.

Figure 4 illustrates the variation in the total number of individuals over time for different values of $R_0$, demonstrating both the disease-free equilibrium (DFE) and endemic equilibrium (EE).

To explore the emergence of limit cycles and periodic behaviors, we conducted further simulations for parameter values leading to oscillatory dynamics. In Figure 5, we examine the variation in the total number of infectious individuals over time for different initial conditions when $R_0>1$, highlighting the stability of endemic equilibrium. This figure shows how the total number of infected individuals changes over time, even under different initial conditions.

Finally, in Figure 6, we present bifurcation diagrams that depict the system’s transitions as parameters such as ${\alpha }_1$, ${\eta }_1$, and ${\sigma }_1$ vary. These diagrams reveal the existence of limit cycles and bifurcation points, providing deeper insights into the dynamics of the infected population.

If ${\eta }_1=0.01234$, ${\alpha }_1=0.02877$, and ${\sigma }_1=0.04000$, then $R_0=0.1837$ (which is less than 1). As a result, the number of susceptible individuals over time initially rises, reaching a peak before leveling off into a steady trend (DFE). On the other hand, the numbers of exposed, unaware infected, aware infected, hospitalized, and recovered individuals all show a declining trend over time, as shown in Figure 4a. If ${\eta }_1=0.1234$, ${\alpha }_1=0.2877$, and ${\sigma }_1=0.4000$, yields $R_0=1.8372$ (which is greater than 1). As a result, the system reaches a state of equilibrium known as endemic equilibrium, shown in Figure 4b. The analysis also explores how parameter variation influences the emergence of limit cycles and complex dynamics. Figure 5a,b show how the fraction of the total number of infected individuals changes over time when $R_0$ is greater than 1. These figures illustrate that for a specific value of $R_0$, the behavior of the number of infected individuals remains consistent, even with different initial numbers of infected individuals. Additionally, the impact of parameters ${\alpha }_1$, ${\eta }_1$, and ${\sigma }_1$ on the infected population ($I_u$) is also depicted in Figure 6, demonstrating how bifurcation curves and limit cycles arise under varying parameter ranges.

11. Conclusions

To summarize, our research provides a novel and enhanced approach to modeling Ebola virus disease (EVD), improving our understanding of disease transmission dynamics and effective countermeasures. The mathematical model introduced in this study incorporates both aware and unaware infected individuals, as well as hospitalized patients, offering a more detailed representation of EVD dynamics. The key findings show that for $R_0<1$, the disease-free equilibrium is stable, meaning the disease will die out, while for $R_0>1$, the endemic equilibrium is stable, indicating that the disease will persist. Sensitivity analysis revealed that the transmission and recovery rates are the most influential factors in controlling $R_0$.

These results have significant implications for public health policy and resource allocation, as they suggest that increasing public awareness and enhancing hospital capacities are critical strategies in managing outbreaks. Additionally, our model’s bifurcation analysis underscores the potential for abrupt changes in disease dynamics based on small parameter variations, further informing intervention strategies. Through numerical simulations, we have demonstrated the practical application of our model, which can guide the optimization of healthcare resources and the development of policies to mitigate the spread of Ebola and similar infectious diseases.