JAM ICCK ICCK Journal of Applied Mathematics 3068-5656 Institute of Central Computation and Knowledge Inc 522 W RIVERSIDE AVE STE N, SPOKANE, WA, 99201, UNITED STATES 10.62762/JAM.2025.550087 Research Article Bifurcation and Stability Analysis of Transmission Dynamics of Ebola Virus Using Seirvh Model https://orcid.org/0000-0002-9657-0328 Ahmad Imtiaz 1 Ullah Wali 1 https://orcid.org/0000-0001-5263-4871 Islam Saeed 2 Ali Nigar 1 Younas Hazrat 1 https://orcid.org/0000-0001-7430-4451 Khan Muhammad Ijaz 2 Department of Mathematics, University of Malakand, Dir Lower, Khyber Pukhtunkhwa, Pakistan Department of Mechanical Engineering, Prince Mohammad Bin Fahd University, Al Khobar 31952, Kingdom of Saudi Arabia Corresponding Author: Saeed Islam. Email: SISLAM@pmu.edu.sa 27 7 2025 1 2 41 51 06 4 2025 16 5 2025 © 2025 by the Authors. Published by Institute of Central Computation and Knowledge. This is an open access article under the CC BY license (https://creativecommons.org/licenses/by/4.0/). 2025 The Authors This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. This article is available from https://www.icck.org/article/abs/jam.2025.550087

This study presents a mathematical framework to analyze the transmission dynamics of the Ebola Virus Disease (EVD) using an extended SEIRVH model. The model incorporates vaccinated and hospitalized compartments, addressing critical factors such as vaccination efficacy, healthcare interventions, and natural disease progression. Differential equations describe the transitions between six population compartments. The study evaluates model stability and bifurcation through well-posedness, positivity, and boundedness analyzes, ensuring realistic and biologically valid solutions. The basic reproduction number, R0, derived from the next generation matrix, serves as a threshold for outbreak control. Local and global stability analyzes of disease-free and endemic equilibria reveal critical insights into epidemic thresholds and long-term dynamics. Furthermore, sensitivity analysis highlights key parameters that influence R0, emphasizing the importance of vaccination and hospitalization in mitigating EVD outbreaks. Numerical simulations validate theoretical findings, underscoring the model's utility in informing effective public health strategies, such as vaccination campaigns and hospitalization measures, for controlling EVD transmission. This research provides a robust analytical and computational tool for understanding and managing the spread of Ebola and similar infectious diseases.

 EBOLA stability analysis sensitivity analysis bifurcation analysis RK-4 method Introduction

Ebola Virus Disease (EVD) is a highly virulent zoonotic illness caused by the Ebola virus, a member of the Filoviridae family. Clinically, EVD presents as an acute hemorrhagic fever syndrome, marked by high-grade fever, severe headache, vomiting, diarrhea, fatigue, and in many cases, internal and external bleeding, which can culminate in multi-organ failure [1, 10]. The disease was first identified in 1976 during simultaneous outbreaks in Sudan and Zaire (now the Democratic Republic of the Congo), with its name derived from the nearby Ebola River [1]. Since then, recurrent outbreaks, especially in Central and West Africa, have underscored its public health importance [2, 4].

The 2014–2016 West African epidemic, which severely impacted Guinea, Liberia, and Sierra Leone, was the largest on record in terms of geographic spread, morbidity, and mortality [4, 25]. EVD outbreaks typically exhibit case-fatality rates ranging from 25% to 90% [4, 10], emphasizing the critical need for timely detection, rapid intervention, and sustainable control strategies.

Transmission occurs primarily via direct contact with bodily fluids—such as blood, feces, or vomit—of symptomatic individuals or deceased patients. Indirect transmission via contaminated fomites has also been documented [3, 10]. In many of the most severely affected regions, weak diagnostic capacity, delayed response, and fragile healthcare infrastructure facilitate rapid amplification of the virus [4, 6].

Mathematical modeling has emerged as an indispensable tool in understanding infectious disease dynamics, particularly for EVD [3, 4, 6, 7, 10]. Compartmental models, especially those based on the classical SEIR (Susceptible–Exposed–Infectious–Recovered) structure, have been widely used to simulate epidemic trajectories, estimate key parameters such as the basic reproduction number R0, and assess the potential impact of public health interventions [16, 17, 21, 22].

In this study, we extend the standard SEIR framework by introducing two additional compartments: vaccinated (V) and hospitalised (H). The V compartment reflects immunity through vaccination, while accounting for imperfect efficacy and waning protection [5, 25]. The H compartment accounts for individuals receiving clinical care in isolation units, recognizing that hospitalisation reduces both mortality and transmission risk [13, 14]. Control measures, such as vaccination, hospitalisation, and awareness-based behavioral change, are incorporated explicitly to evaluate their respective effects on epidemic outcomes [11, 12, 14, 15].

The resulting model is a nonlinear system of ordinary differential equations (ODEs) with bilinear incidence, designed to capture both biological and policy-driven transitions. Below we describe the model formulation.

 Model Formulation

To explore the transmission and control dynamics of EVD, we partition the total human population into six mutually exclusive epidemiological compartments:

S⁢(t): Susceptible individuals.

E⁢(t): Exposed (infected but not yet infectious) individuals.

I⁢(t): Infectious individuals.

V⁢(t): Vaccinated individuals.

H⁢(t): Hospitalized individuals under clinical care.

R⁢(t): Recovered individuals with acquired immunity.

Model Assumptions

The model is governed by the following assumptions, consistent with established EVD dynamics [3, 4, 5, 9]:

Recruitment into the susceptible population occurs at a constant rate Ah.

Susceptible individuals are vaccinated at rate α1; vaccine-induced immunity wanes at rate α2.

Susceptible individuals contract EVD via effective contact with infectious individuals at rate β, modeled with a bilinear incidence term β⁢S⁢I [3, 4].

Exposed individuals progress to the infectious class at rate δ1, or are vaccinated (e.g., ring vaccination or post-exposure prophylaxis) at rate δ2 [5, 25].

Vaccinated individuals may still become infected due to vaccine failure at rate δ3.

Infectious individuals recover at rate γ1 or are hospitalized at rate γ2 [13].

Hospitalized individuals recover or die (exit the compartment) at rate γ3.

All compartments are subject to a uniform natural mortality or removal rate k.

Descriptions of parameters used in the model.

Parameter Interpretation
A h Recruitment rate into the susceptible class.
α 1 Vaccination rate for susceptible individuals.
α 2 Rate of waning immunity among vaccinated individuals.
β Effective contact rate between susceptible and infectious individuals.
δ 1 Progression rate from exposed to infectious.
δ 2 Vaccination/intervention rate among exposed individuals.
δ 3 Vaccine failure rate.
γ 1 Recovery rate of infectious individuals.
γ 2 Hospitalisation rate of infectious individuals.
γ 3 Recovery/mortality rate in hospitalized individuals.
k Natural mortality/removal rate across all compartments.

Mathematical Model

The dynamics of the EVD transmission system are captured by the following set of nonlinear ODEs, with parameters defined in Table 1:

{ d ⁢ S d ⁢ t = A h + α 2 ⁢ V − ( α 1 + k + β ⁢ I ) ⁢ S , d ⁢ E d ⁢ t = β ⁢ S ⁢ I − ( δ 1 + δ 2 + k ) ⁢ E , d ⁢ I d ⁢ t = δ 1 ⁢ E + δ 3 ⁢ V − ( γ 1 + γ 2 + k ) ⁢ I , d ⁢ V d ⁢ t = α 1 ⁢ S + δ 2 ⁢ E − ( δ 3 + α 2 + k ) ⁢ V , d ⁢ H d ⁢ t = γ 2 ⁢ I − ( γ 3 + k ) ⁢ H , d ⁢ R d ⁢ t = γ 1 ⁢ I + γ 3 ⁢ H − k ⁢ R .

Subject to the following initial conditions:

S ⁢ ( 0 ) > 0 , E ⁢ ( 0 ) ≥ 0 , I ⁢ ( 0 ) ≥ 0 , V ⁢ ( 0 ) ≥ 0 , H ⁢ ( 0 ) ≥ 0 , R ⁢ ( 0 ) ≥ 0 .

Figure 1 presents the flow diagram representing the compartmental structure and transitions of the proposed model. Each parameter and transition is biologically motivated and derived from literature [3, 4, 10].

Flow diagram of the proposed Ebola virus transmission model.

 Theoretical Analysis Positivity of Solutions <bold>.</bold>

Given non-negative initial conditions, all solutions of the system (1) remain non-negative for all t≥0.

Proof..

We prove the non-negativity of each compartment in system (1) using contradiction and differential inequalities.

Step 1: S⁢(t)≥0Suppose there exists a first time t0>0 such that S⁢(t0)=0, with S⁢(t)>0 for t<t0, and S⁢(t)<0 for some t>t0. The governing equation is:

d ⁢ S d ⁢ t = A h + α 2 ⁢ V − ( α 1 + k + β ⁢ I ) ⁢ S .

At t=t0, since S⁢(t0)=0, we have:

d ⁢ S d ⁢ t | t = t 0 = A h + α 2 ⁢ V ⁢ ( t 0 ) ≥ 0 ,

which contradicts S⁢(t)<0 for t>t0. Thus, S⁢(t)≥0 for all t≥0.

Step 2: E⁢(t)≥0

d ⁢ E d ⁢ t = β ⁢ S ⁢ I − ( δ 1 + δ 2 + k ) ⁢ E .

If E⁢(t0)=0, then d⁢Ed⁢t=β⁢S⁢I≥0, so E⁢(t)≥0.

Step 3: I⁢(t)≥0

d ⁢ I d ⁢ t = δ 1 ⁢ E + δ 3 ⁢ V − ( γ 1 + γ 2 + k ) ⁢ I .

If I⁢(t0)=0, then d⁢Id⁢t=δ1⁢E+δ3⁢V≥0.

Step 4: V⁢(t)≥0

d ⁢ V d ⁢ t = α 1 ⁢ S + δ 2 ⁢ E − ( δ 3 + α 2 + k ) ⁢ V .

If V⁢(t0)=0, then d⁢Vd⁢t=α1⁢S+δ2⁢E≥0.

Step 5: H⁢(t)≥0

d ⁢ H d ⁢ t = γ 2 ⁢ I − ( γ 3 + k ) ⁢ H .

If H⁢(t0)=0, then d⁢Hd⁢t=γ2⁢I≥0.

Step 6: R⁢(t)≥0

d ⁢ R d ⁢ t = γ 1 ⁢ I + γ 3 ⁢ H − k ⁢ R .

If R⁢(t0)=0, then d⁢Rd⁢t=γ1⁢I+γ3⁢H≥0.

Conclusion: All compartments remain non-negative for all t≥0. ∎

 Disease-Free Equilibrium Point

The disease-free equilibrium (DFE) represents a state where no infection exists in the population:

( S ∗ , E ∗ , I ∗ , V ∗ , H ∗ , R ∗ ) = ( A h α 1 + k , 0 , 0 , 0 , 0 , 0 ) .

The DFE is stable if the basic reproduction number R0≤1. When R0>1, the disease may invade and persist in the population.

 Endemic Equilibrium Point

An endemic equilibrium refers to a state where the disease persists over time. The number of infections, recoveries, and deaths balance such that:

S ∗ = A h − V ∗ α 1 + k + β ⁢ I ∗ , E ∗ = β ⁢ S ∗ ⁢ I ∗ δ 1 + δ 2 + k , I ∗ = δ 1 ⁢ E ∗ + δ 3 ⁢ V ∗ γ 1 + γ 2 + k , V ∗ = α 1 ⁢ S ∗ + δ 2 ⁢ E ∗ δ 3 + α 2 + k , H ∗ = γ 2 ⁢ I ∗ γ 3 + k , R ∗ = γ 1 ⁢ I ∗ + γ 3 ⁢ H ∗ k .

Model Analysis and Boundedness

Define the total human population as:

N ⁢ ( t ) = S ⁢ ( t ) + E ⁢ ( t ) + I ⁢ ( t ) + V ⁢ ( t ) + H ⁢ ( t ) + R ⁢ ( t ) .

Differentiating, we obtain:

d ⁢ N d ⁢ t = A h − k ⁢ N ⁢ ( t ) .

Solving the differential equation yields:

N ⁢ ( t ) = N ⁢ ( 0 ) ⁢ e − k ⁢ t + A h k ⁢ ( 1 − e − k ⁢ t ) .

Hence, as t→∞, the population tends to:

lim t → ∞ N ⁢ ( t ) = A h k .

<bold>.</bold>

The region

Ω = { ( S , E , I , V , H , R ) ∈ ℝ + 6 : N ⁢ ( t ) ≤ A h k }

is positively invariant.

Proof..

From the total population dynamics:

d ⁢ N d ⁢ t = A h − k ⁢ N ,

with the integrating factor method, we obtain:

N ⁢ ( t ) ≥ N ⁢ ( 0 ) ⁢ e − k ⁢ t + A h k ⁢ ( 1 − e − k ⁢ t ) .

Therefore, N⁢(t)→Ahk as t→∞, and solutions are bounded within region Ω. Furthermore, since each component satisfies:

( S ⁢ ( t ) E ⁢ ( t ) I ⁢ ( t ) V ⁢ ( t ) H ⁢ ( t ) R ⁢ ( t ) ) ∈ ℝ + 6 ,

the model is mathematically well-posed and epidemiologically meaningful. ∎

 Basic Reproduction Number (<inline-formula><mml:math alttext="R_{0}" display="inline"><mml:msub><mml:mi>R</mml:mi><mml:mn>0</mml:mn></mml:msub></mml:math></inline-formula>)

The basic reproduction number, denoted as R0, is a critical threshold used to evaluate the potential for disease spread in a fully susceptible population. Specifically, R0 represents the expected number of secondary infections produced by a single infectious individual introduced into a completely susceptible population. When R0<1, the disease cannot invade the population and will eventually die out. Conversely, if R0>1, the disease may cause an outbreak and persist in the population. Thus, understanding and computing R0 is crucial for predicting the future dynamics of the disease and guiding effective control strategies.

Numerous studies have addressed the determination of R0 for various epidemiological models [15, 16]. In this work, we compute R0 for the proposed model (1) using the next-generation matrix method [17, 18]. The infection subsystem is isolated, and the matrices F and V are constructed in accordance with the formulation given in [8].

The infection subsystem of the model (1) is given by:

{ d ⁢ E d ⁢ t = β ⁢ S ⁢ I − ( δ 1 + δ 2 + k ) ⁢ E , d ⁢ I d ⁢ t = δ 1 ⁢ E + δ 3 ⁢ V − ( γ 1 + γ 2 + k ) ⁢ I , d ⁢ V d ⁢ t = α 1 ⁢ S + δ 2 ⁢ E − ( δ 3 + α 2 + k ) ⁢ V .

The new infection matrix F and its Jacobian at the disease-free equilibrium are:

F = ( β ⁢ S ⁢ I 0 0 ) , F ∗ = ( 0 β ⁢ S 0 0 0 0 0 0 0 0 ) .

The transition matrix V is defined as:

V = ( ( δ 1 + δ 2 + k ) ⁢ E ( γ 1 + γ 2 + k ) ⁢ I − δ 1 ⁢ E − δ 3 ⁢ V ( δ 3 + α 2 + k ) ⁢ V − δ 2 ⁢ E ) .

The inverse of the Jacobian of V evaluated at the disease-free equilibrium is:

V ∗ − 1 = ( δ 1 + δ 2 + k 0 0 − δ 1 γ 1 + γ 2 + k − δ 3 − δ 2 0 δ 3 + α 2 + k ) .

Multiplying F∗ and V∗−1 yields:

F ⁢ V − 1 = ( β ⁢ S 0 ⁢ δ 1 ( δ 1 + δ 2 + k ) ⁢ ( γ 1 + γ 2 + k ) β ⁢ S 0 γ 1 + γ 2 + k 0 0 0 0 0 0 0 ) .

The basic reproduction number R0 is given by the spectral radius of F⁢V−1:

R 0 = β ⁢ δ 1 ⁢ A h ( δ 1 + δ 2 + k ) ⁢ ( γ 1 + γ 2 + k ) ⁢ ( α 1 + k ) .

The effects of the parameters β, Ah, and k on R0 are illustrated in Figure 2, highlighting their influence on the potential for disease spread.

Effects of β, Ah, and k on R0.

Local Stability of the Disease-Free Equilibrium <bold>.</bold>

The disease-free equilibrium E0 of system (1) is locally asymptotically stable if R0<1 and unstable if R0>1.

Proof..

At the disease-free equilibrium point E0, the Jacobian matrix of system (1) is:

J ⁢ E 0 = ( − ( α 1 + k ) 0 β ⁢ S ∗ α 2 0 0 0 − Z 1 β ⁢ S ∗ 0 0 0 0 δ 1 − Z 2 δ 3 0 0 α 1 δ 2 0 − ( δ 3 + α 2 + k ) 0 0 0 0 γ 2 0 − ( γ 3 + k ) 0 0 0 γ 1 0 γ 3 − k ) ,

where Z1=δ1+δ2+k and Z2=γ1+γ2+k. The eigenvalues of this matrix are:

λ 1 = − ( α 1 + k ) , λ 2 = − ( γ 1 + γ 2 + k ) , λ 3 = − k , λ 4 = − Z 1 , λ 5 = − Z 1 ⁢ Z 2 , λ 6 = − Z 1 ⁢ Z 2 ⁢ ( δ 3 + α 2 + k ) .

Since all eigenvalues are negative when R0<1, the disease-free equilibrium is locally asymptotically stable by the Routh-Hurwitz criterion [21, 22, 24]. ∎

 Global Stability of the Endemic Equilibrium <bold>.</bold>

If R0>1, then the endemic equilibrium E1 of system (1) is globally asymptotically stable in Ω.

Proof..

For R0>1, the existence of a unique endemic equilibrium E1 is guaranteed. We consider the standard Lyapunov function [23]:

V ⁢ ( x ) = ∑ i = 1 n c i 2 ⁢ ( x i − x i ∗ ) 2 ,

which, for our system, takes the form:

V ( S , E , I , V , H , R ) = 1 2 [ ( S − S ∗ ) + ( E − E ∗ ) + ( I − I ∗ ) + ( V − V ∗ ) + ( H − H ∗ ) + ( R − R ∗ ) ] 2 .

Differentiating (10) along solutions of system (1), and noting that:

d d ⁢ t ⁢ ( S + E + I + V + H + R ) = − k ⁢ ( S + E + I + V + H + R ) ,

we obtain:

d ⁢ V d ⁢ t = − k ( S + E + I + V + H + R ) ⋅ [ ( S − S ∗ ) + ( E − E ∗ ) + ( I − I ∗ ) + ( V − V ∗ ) + ( H − H ∗ ) + ( R − R ∗ ) ] .

Clearly, d⁢Vd⁢t≤0, and equality holds if and only if the state variables equal their endemic equilibrium values. Hence, by LaSalle's Invariance Principle [23], the endemic equilibrium E1 is globally asymptotically stable in Ω. ∎

 Sensitivity Analysis of <inline-formula><mml:math alttext="\bm{R_{0}}" display="inline"><mml:msub><mml:mi>𝑹</mml:mi><mml:mn>𝟎</mml:mn></mml:msub></mml:math></inline-formula>

The sensitivity analysis of the basic reproduction number R0 with respect to model parameters plays a crucial role in understanding the influence of each parameter on disease transmission, control, and treatment effectiveness, as demonstrated in similar infectious disease models [19]. This analysis identifies which parameters most significantly impact R0 and thus should be prioritized in disease mitigation strategies.

We examine the sensitivity of R0 in relation to the key parameters of the Ebola virus model (1). The sensitivity index of a parameter ψ is computed using the normalized forward sensitivity index as proposed by Chitnis et al. [21], defined as:

X ψ R 0 = ∂ R 0 ∂ ψ × ψ R 0

The sensitivity indices for relevant model parameters are computed as follows, with their impacts illustrated in Figure 3:

X β R 0 = 1 , X A h R 0 = 1 , X δ 1 R 0 = − δ 1 δ 1 + δ 2 + k , X δ 2 R 0 = − δ 2 δ 1 + δ 2 + k , X k R 0 = − k ⁢ ( 1 δ 1 + δ 2 + k + 1 γ 1 + γ 2 + k + 1 α 1 + k ) , X γ 1 R 0 = − γ 1 ⁢ ( ( δ 1 + δ 2 + γ 1 ) ⁢ ( α 1 + k ) + ( γ 1 + γ 2 + k ) ⁢ ( α 1 + k ) ( δ 1 + δ 2 + γ 1 ) ⁢ ( γ 1 + γ 2 + k ) ⁢ ( α 1 + k ) ) , X γ 2 R 0 = − γ 2 ⁢ ( ( δ 1 + δ 2 + γ 1 ) ⁢ ( α 1 + k ) + ( γ 1 + γ 2 + k ) ⁢ ( δ 1 + δ 2 + γ 1 ) ( δ 1 + δ 2 + γ 1 ) ⁢ ( γ 1 + γ 2 + k ) ⁢ ( α 1 + k ) ) , X α 1 R 0 = − α 1 ⋅ β ⁢ A h ⁢ ( δ 1 + δ 2 + γ 1 ) ( δ 1 + δ 2 + γ 1 ) 2 ⁢ ( γ 1 + γ 2 + k ) 2

Sensitivity indices of parameters influencing R0.

 Optimal Control System

Dynamics of EVD compartments with and without control interventions (m1, m2).

To further mitigate the Ebola Virus Disease (EVD), we expand the model (1) by incorporating two time-dependent control functions: m1⁢(t) and m2⁢(t). Here, m1⁢(t) represents enhanced diagnostic efforts for exposed individuals, while m2⁢(t) captures the immediate treatment of infected individuals. The extended model is given by:

d ⁢ S d ⁢ t = A h + α 2 ⁢ V − ( α 1 + k ) ⁢ S − ( 1 − m 1 ) ⁢ β ⁢ I ⁢ S , d ⁢ E d ⁢ t = ( 1 − m 1 ) ⁢ β ⁢ I ⁢ S − ( δ 1 + δ 2 + k ) ⁢ E , d ⁢ I d ⁢ t = δ 1 ⁢ E + δ 3 ⁢ V − ( γ 1 + m 2 ⁢ γ 2 + k ) ⁢ I , d ⁢ V d ⁢ t = α 1 ⁢ S + δ 2 ⁢ E − ( δ 3 + α 2 + k ) ⁢ V , d ⁢ H d ⁢ t = m 2 ⁢ γ 2 ⁢ I − ( γ 3 + k ) ⁢ H , d ⁢ R d ⁢ t = γ 1 ⁢ I + γ 3 ⁢ H − k ⁢ R .

The impact of these control measures on the dynamics of the EVD compartments is illustrated in Figure 4. The goal is to minimize the number of exposed and infected individuals while also minimizing the costs of the control measures. This is achieved via the following cost functional:

𝒥 ⁢ ( m 1 , m 2 ) = ∫ t 0 t f [ ℬ 1 ⁢ E + ℬ 2 ⁢ I ] ⁢ 𝑑 t + 1 2 ⁢ ( ℬ 3 ⁢ m 1 2 ⁢ ( t ) + ℬ 4 ⁢ m 2 2 ⁢ ( t ) ) ,

subject to the admissible control set:

𝒰 = { m j ⁢ ( t ) : 0 ≤ m j ⁢ ( t ) ≤ 1 , j = 1 , 2 , t ∈ [ 0 , t f ] } .

Characterization of Optimal Control

Applying Pontryagin's Maximum Principle [23], a widely used approach in infectious disease control models [20], the Hamiltonian ℋ for the system is given as:

ℋ = ℬ 1 ⁢ E + ℬ 2 ⁢ I + 1 2 ⁢ ( ℬ 3 ⁢ m 1 2 + ℬ 4 ⁢ m 2 2 ) + ∑ i = 1 6 λ i ⁢ d ⁢ x i d ⁢ t ,

where λi are the adjoint variables corresponding to state variables S,E,I,V,H, and R.

The adjoint system is governed by:

d ⁢ λ 1 d ⁢ t = λ 1 ⁢ ( α 1 + k ) + ( λ 1 − λ 2 ) ⁢ ( 1 − m 1 ) ⁢ β ⁢ I + λ 4 ⁢ α 1 , d ⁢ λ 2 d ⁢ t = − ℬ 1 + ( δ 1 + δ 2 + k ) ⁢ λ 2 − λ 3 ⁢ δ 1 − λ 4 ⁢ δ 2 , d ⁢ λ 3 d ⁢ t = − ℬ 2 + ( λ 1 − λ 2 ) ⁢ ( 1 − m 1 ) ⁢ β ⁢ S + λ 3 ⁢ ( γ 1 + m 2 ⁢ γ 2 + k ) − λ 5 ⁢ m 2 ⁢ γ 2 − λ 6 ⁢ γ 1 , d ⁢ λ 4 d ⁢ t = − λ 1 ⁢ α 1 − λ 3 ⁢ δ 3 + λ 4 ⁢ ( δ 2 + α 2 + k ) , d ⁢ λ 5 d ⁢ t = λ 5 ⁢ ( γ 3 + k ) − λ 6 ⁢ γ 3 , d ⁢ λ 6 d ⁢ t = λ 6 ⁢ k ,

with transversality conditions:

λ i ⁢ ( t f ) = 0 , for ⁢ i = 1 , … , 6 .

The optimal controls are characterized by:

m 1 ∗ ⁢ ( t ) = min ⁡ ( 1 , max ⁡ ( 0 , ( λ 2 − λ 1 ) ⁢ β ⁢ I ⁢ S ℬ 3 ) ) , m 2 ∗ ⁢ ( t ) = min ⁡ ( 1 , max ⁡ ( 0 , ( λ 3 − λ 5 ) ⁢ γ 2 ⁢ I ℬ 4 ) ) .

Combining the state system (1), the adjoint system (4), and the optimal controls (5), we obtain the full optimality system for simulation and control implementation.

 Discussion

The relationship between the infection rate β and the basic reproduction number R0 is direct and proportional. An increase in β leads to a corresponding increase in R0, indicating that transmission intensity escalates with higher infection rates. Therefore, reducing β is essential for outbreak control. This can be achieved through interventions such as immunization, isolation, public health education, and hygiene improvements.

The sensitivity index of δ1 is negative, implying that an increase in δ1 reduces R0. A higher mortality rate in the first infectious class shortens the infectious period, thereby limiting disease transmission. Similarly, an increase in δ2 also lowers R0, underscoring the role of mortality in reducing the number of secondary infections.

For the recovery rate k, its sensitivity index is also negative. Enhancing the recovery rate through medical treatment or supportive care diminishes the number of infectious individuals, hence reducing R0.

The transition parameters γ1 and γ2 represent movement between different risk groups or health statuses. Their negative sensitivity indices suggest that facilitating transitions (e.g., hospitalization or effective isolation) can mitigate disease spread.

Furthermore, the natural death rate α1 exhibits a negative sensitivity index. A higher α1 reduces the pool of susceptible individuals, thereby lowering the overall transmission potential.

 Numerical Simulation

To analyze the impact of model parameters on the disease dynamics, we evaluate the sensitivity of the basic reproduction number R0 using the normalized forward sensitivity index defined by Chitnis et al. [21], given by Equation (13).

Γ θ R 0 = ∂ R 0 ∂ θ ⋅ θ R 0 ,

where θ denotes a given model parameter. Sensitivity analysis helps identify key parameters for intervention strategies.

Dynamics of compartments S,E,I,V,H, and R.

Effect on Susceptible class.

Effect of hospitalization on susceptible population.

Effect on exposed class.

Effect on infectious class.

Effect on recovered class.

Vaccination impact on susceptible class.

Vaccination impact on exposed class.

Vaccination impact on infected class.

Vaccination impact on recovered class.

Hospitalized population under vaccination.

Figure 5 illustrates the time evolution of each compartment in the SEIRVH model. During simulations, a subset of parameters is varied while others remain constant to assess system response. When β=0.003, we find that R0=0.9<1, and the disease-free equilibrium (DFE) is locally asymptotically stable.

Similarly, for β=0.0006, R0=0.1499<1, the system also converges to the DFE. In contrast, when β=0.006, resulting in R0=1.4990>1, the system approaches an endemic equilibrium, aligning with Theorem 4.1 and Theorem 4.2.

Figures 6–15 show that increasing β accelerates the depletion of the susceptible population and intensifies the epidemic. A smaller β (e.g., 0.1) leads to slower disease spread, while a higher β (e.g., 0.5) results in rapid increases in exposed and infectious individuals.

Vaccination effects are evident in Figures 11–14. With vaccination, the susceptible, exposed, and infectious populations decline more quickly. The number of recovered individuals increases, indicating improved control over disease spread.

 Conclusion Remarks

This study presents a comprehensive analytical and numerical investigation of the Ebola Virus Disease (EVD) through an extended SEIRVH model that incorporates vaccination and hospitalization. The model's well-posedness, positivity, and boundedness affirm its biological validity.

Using the next-generation matrix approach, we derived the basic reproduction number R0 and assessed local and global stability for both disease-free and endemic equilibria. The Lyapunov function method confirms the global stability conditions, while sensitivity analysis reveals the most influential parameters on R0.

A backward bifurcation was identified using Center Manifold theory, emphasizing the challenge of achieving disease eradication even when R0<1. Optimal control analysis incorporating vaccination, therapy, and public awareness further illustrates effective mitigation strategies.

Numerical simulations validate the theoretical results and demonstrate the crucial role of reducing β through public health interventions. Moreover, the Non-standard Finite Difference (NSFD) scheme outperformed the classical RK-4 method in simulating disease dynamics.

Overall, this model provides a robust tool for understanding and managing EVD outbreaks. It highlights the necessity of integrated strategies combining vaccination, hospitalization, and timely interventions, offering practical guidance for public health authorities to mitigate the devastating impacts of infectious diseases.

 Acknowledgments

This work was supported without any funding.

 Conflict of interest

The authors declare no conflicts of interest.

 References Piot, P., Sureau, P., Breman, G., Heymann, D., Kintoki, V., Masamba, M., & Mbuy, M. (1978). Clinical aspects of Ebola virus infection in Yambuku area, Zaire 1976. African Journal of Infectious Diseases, 210(96), 34–51. Chretien, J. P., Riley, S., & George, D. B. (2015). Mathematical modeling of the West Africa Ebola epidemic. Elife, 4, e09186. [https://doi.org/10.7554/eLife.09186] Durojaye, M., & Ajie, J. (2017). Mathematical model of the spread and control of Ebola virus disease. Applied Mathematics, 7, 23–31. Mbah, G. C., Sunday, O. I., Ojoma, A. Q., Christopher, A., & Chukwudi, O. (2023). Mathematical modelling approach of the study of Ebola virus disease transmission dynamics in a developing country. African Journal of Infectious Diseases, 17(1), 10–26. Tadmon, C., & Kengne, J. N. (2022). Mathematical analysis of a model of Ebola disease with control measures. International Journal of Biomathematics, 15(07), 2250048. [https://doi.org/10.1142/S1793524522500486] Djiomba Njankou, S. D. (2019). Mathematical models of Ebola virus disease with socio-economic dynamics (Doctoral dissertation, Stellenbosch: Stellenbosch University). Ahmad, W., & Abbas, M. (2021). Effect of quarantine on transmission dynamics of Ebola virus epidemic: a mathematical analysis. The European Physical Journal Plus, 136(4), 1-33. [https://doi.org/10.1140/epjp/s13360-021-01360-9] La Salle, J. P. (1976). The stability of dynamical systems. Society for Industrial and Applied Mathematics. Shah, I., Alrabaiah, H., & Ozdemir, B. (2023). Using advanced analysis together with fractional order derivative to investigate a smoking tobacco cancer model. Results in Physics, 51, 106700. [https://doi.org/10.1016/j.rinp.2023.106700] Juga, M. L., Nyabadza, F., & Chirove, F. (2021). An Ebola virus disease model with fear and environmental transmission dynamics. Infectious Disease Modelling, 6, 545-559. [https://doi.org/10.1016/j.idm.2021.03.002] Funk, S., Gilad, E., Watkins, C., & Jansen, V. A. A. (2009). The spread of awareness and its impact on epidemic outbreaks. Proceedings of the National Academy of Sciences, 106(16), 6872–6877. [https://doi.org/10.1073/pnas.0810762106] Chamchod, F., & Britton, N. F. (2012). On the dynamics of a two-strain influenza model with isolation. Mathematical Modelling of Natural Phenomena, 7(3), 49–61. [https://doi.org/10.1051/mmnp/20127305] Day, T., Park, A., Madras, N., Gumel, A. B., & Wu, J. (2006). When is quarantine a useful control strategy for emerging infectious diseases? American Journal of Epidemiology, 163(5), 479–485. [https://doi.org/10.1093/aje/kwj056] Safi, M., & Gumel, A. B. (2013). Dynamics of a model with quarantine-adjusted incidence and quarantine of susceptible individuals. Journal of Mathematical Analysis and Applications, 399(2), 565–575. [https://doi.org/10.1016/j.jmaa.2012.10.015] Zakary, O., Rachik, M., & Elmouki, I. (2017). On the analysis of a multi-regions discrete SIR epidemic model: An optimal control approach. International Journal of Dynamics and Control, 5(3), 917–930. [https://doi.org/10.1007/s40435-016-0233-2] Brauer, F., Castillo-Chavez, C., & Castillo-Chavez, C. (2012). Mathematical models in population biology and epidemiology (Vol. 2, No. 10). New York: springer. Chavez, C. C., Feng, Z., & Huang, W. (2002). On the computation of R0 and its role on global stability. Mathematical approaches for emerging and re-emerging infection diseases: An introduction, 125, 31-65. Thabet, S. T., Abdo, M. S., Shah, K., & Abdeljawad, T. (2020). Study of transmission dynamics of COVID-19 mathematical model under ABC fractional order derivative. Results in Physics, 19, 103507. [https://doi.org/10.1016/j.rinp.2020.103507] Aja, R. O., Omale, D., & Mbah, G. C. E. (2017). Sensitivity Analysis of the Mathematical Model on the Control of HBV-HDV co-infection Transmission Dynamics in a Given Population. Journal of the Nigerian Association of Mathematical Physics, 39, 457-470. Kotola, B. S., Teklu, S., & Abebaw, Y. F. (2023). Bifurcation and optimal control analysis of HIV/AIDS and COVID-19 co-infection model with numerical simulation. PLoS ONE, 18(5), e0284759. [https://doi.org/10.1371/journal.pone.0284759] Chitnis, N., Hyman, J. M., & Cushing, J. M. (2008). Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model. Bulletin of Mathematical Biology, 70(5), 1272–1296. [https://doi.org/10.1007/s11538-008-9299-0] Van den Driessche, P., & Watmough, J. (2008). Further notes on the basic reproduction number. Journal of Mathematical Biology, 56(1), 159–178. [https://doi.org/10.1007/978-3-540-78911-6_6] Delamater, P. L., Street, E. J., Leslie, T. F., Yang, Y. T., & Jacobsen, K. H. (2019). Complexity of the basic reproduction number R0. Emerging Infectious Diseases, 25(1), 1–4. [https://doi.org/10.3201/eid2501.171901] Mahardika, R., Widowati, & Sumanto, Y. D. (2019, May). Routh-hurwitz criterion and bifurcation method for stability analysis of tuberculosis transmission model. In Journal of physics: Conference series (Vol. 1217, No. 1, p. 012056). IOP Publishing. [https://doi.org/10.1088/1742-6596/1217/1/012056] Brettin, A., Goldthorpe, R., Weishaar, K., & Erovenko, I. V. (2018). Ebola could be eradicated through voluntary vaccination. Royal Society Open Science, 5(5), 171591. [https://doi.org/10.1098/rsos.171591]