Modeling and Analysis of SIRR Model (Ebola Transmission Dynamics Model) with Delay Differential Equation

Akinleye Emmanuel Lasekan; Joshua Oluwasegun Agbomola; Kabir Oluwatobi Idowu; Babatunde Ademola Kannike; Esther Oluwatoyin Mulero; Temitope Senami Gandonu; Solari Myrjuari Elee

doi:10.12688/f1000research.168361.1

Home Browse Modeling and Analysis of SIRR Model (Ebola Transmission Dynamics Model)...

ALL Metrics

-

Views

-

Downloads

Get PDF

Get XML

Export

▬

✚

Research Article

Modeling and Analysis of SIRR Model (Ebola Transmission Dynamics Model) with Delay Differential Equation

[version 1; peer review: 2 approved]

Akinleye Emmanuel Lasekan¹, Joshua Oluwasegun Agbomola², Kabir Oluwatobi Idowu ³, [...] Babatunde Ademola Kannike¹, Esther Oluwatoyin Mulero¹, Temitope Senami Gandonu², Solari Myrjuari Elee⁴

Akinleye Emmanuel Lasekan¹, Joshua Oluwasegun Agbomola², [...] Kabir Oluwatobi Idowu ³, Babatunde Ademola Kannike¹, Esther Oluwatoyin Mulero¹, Temitope Senami Gandonu², Solari Myrjuari Elee⁴

PUBLISHED 02 Sep 2025

Author details Author details

¹ Department of Mathematics, Lagos State University, Ojo, Lagos, Nigeria
² Department of Mathematics, Tulane University, New Orleans, Louisiana, USA
³ Department of Mathematics, Purdue University, West Lafayette, Indiana, USA
⁴ Department of Mathematics, University of Delaware, Newark, Delaware, USA

Akinleye Emmanuel Lasekan
Roles: Conceptualization, Formal Analysis, Writing – Original Draft Preparation

Joshua Oluwasegun Agbomola
Roles: Methodology, Validation, Writing – Original Draft Preparation

Kabir Oluwatobi Idowu
Roles: Conceptualization, Formal Analysis, Methodology, Supervision, Writing – Review & Editing

Babatunde Ademola Kannike
Roles: Project Administration, Writing – Original Draft Preparation, Writing – Review & Editing

Esther Oluwatoyin Mulero
Roles: Writing – Original Draft Preparation, Writing – Review & Editing

Temitope Senami Gandonu
Roles: Writing – Original Draft Preparation, Writing – Review & Editing

Solari Myrjuari Elee
Roles: Writing – Original Draft Preparation, Writing – Review & Editing

OPEN PEER REVIEW

REVIEWER STATUS

This article is included in the Pathogens gateway.

This article is included in the Emerging Diseases and Outbreaks gateway.

Abstract

Background

Ebola virus disease (EVD) is a severe and often fatal illness with high transmission potential and recurring outbreaks. Traditional compartmental models often neglect biologically important delays, such as the latent period before an infected individual becomes infectious, limiting their ability to capture real-world epidemic patterns. Including such delays can provide a more accurate understanding of outbreak persistence and control strategies.

Methods

In this study, we develop and analyze a novel deterministic SIRR model that captures the complex transmission dynamics of Ebola by explicitly combining nonlinear incidence rates with a delay differential equation framework. Unlike traditional models, this approach integrates a biologically motivated delay to represent the latent period before infectiousness, providing a more realistic depiction of disease spread. The basic reproduction number (R₀) is derived using the next-generation matrix, and local stability for disease-free and endemic equilibria is established. Using center manifold theory, we investigate transcritical bifurcation at R₀ = 1, while Hopf bifurcation analysis determines when delays trigger oscillatory epidemics. Sensitivity analysis identifies parameters most influencing R₀, and numerical simulations are performed using the fourth-order Runge–Kutta method.

Results

The main novelty of this work lies in its detailed investigation of how delays influence outbreak persistence and can trigger oscillatory epidemics, patterns often observed in practice but rarely captured by classic models. For R₀< 1, the disease-free equilibrium is locally asymptotically stable; for R₀> 1, an endemic equilibrium emerges. Increasing delays destabilizes the system, amplifying peak infections, prolonging outbreaks, and producing sustained oscillations. Isolation of recovered individuals (c) significantly reduces R_0, while transmission rate (β), recruitment rate (Λ), and isolation transition rate (ρ) are identified as the most sensitive parameters.

Conclusions

Accounting for delayed recovery dynamics is crucial for accurately predicting outbreak patterns and designing effective interventions. This delay-based, nonlinear-incidence model offers a robust analytical and computational framework for guiding public health strategies, with direct implications for reducing transmission, shortening outbreak duration, and preventing epidemic resurgence.

Keywords

SIRR model, delay differential equation, Equilibrium point \sep stability, center manifold theory, Bifurcation Analysis

Corresponding author: Kabir Oluwatobi Idowu

Competing interests: No competing interests were disclosed.

Grant information: The author(s) declared that no grants were involved in supporting this work.

Copyright: © 2025 Lasekan AE et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

How to cite: Lasekan AE, Oluwasegun Agbomola J, Idowu KO et al. Modeling and Analysis of SIRR Model (Ebola Transmission Dynamics Model) with Delay Differential Equation [version 1; peer review: 2 approved]. F1000Research 2025, 14:857 (https://doi.org/10.12688/f1000research.168361.1) First published: 02 Sep 2025, 14:857 (https://doi.org/10.12688/f1000research.168361.1) Latest published: 02 Sep 2025, 14:857 (https://doi.org/10.12688/f1000research.168361.1)

1. Introduction

Epidemiological modeling plays a pivotal role in understanding and controlling the spread of infectious diseases.¹ These models enable researchers and policymakers to analyze disease transmission dynamics, predict outbreaks, and design effective intervention strategies. Among the earliest and most influential frameworks is the SIR model, introduced by Kermack and McKendrick in 1927.¹ This simple yet profound model divided the population into three compartments: susceptible, infected, and recovered, and demonstrated the concept of a threshold level of population immunity required to stop an epidemic. The SIR model not only explained how epidemics unfold but also underscored the importance of early interventions. Over the decades, it has remained a cornerstone of infectious disease modeling and public health strategies, forming the basis for countless extensions and refinements in mathematical epidemiology.^1,2

In the mid to late 20th century, the field underwent significant expansion, driven by the work of researchers like Anderson and May.³ They broadened the scope of epidemiological models by incorporating critical factors such as heterogeneous mixing patterns, disease vectors, and variable transmission rates. These innovations made the models more applicable to real-world scenarios, providing insights into the dynamics of diseases like malaria, HIV, and schistosomiasis. Their studies emphasized the importance of understanding transmission dynamics, host-vector relationships, and the role of immunity in shaping the persistence and control of endemic diseases. Additionally, they highlighted the potential impact of control measures such as vaccination and vector management on disease eradication efforts.⁴

Building on these advancements, Hethcote⁵ made substantial contributions by addressing nonlinear dynamics, seasonality, and spatial heterogeneity in epidemiological models. His work introduced new dimensions to understanding disease outbreaks, particularly in cases where environmental or seasonal factors play a significant role. These refinements underscored the complexity of real-world disease transmission and the need for models capable of capturing multifaceted interactions among biological, social, and environmental variables. The concept of the basic reproduction number R₀, as discussed by Heesterbeek and Dietz,⁶ further formalized how to quantify the potential of a pathogen to spread in a population.

The urgency and complexity of Ebola virus outbreaks have driven widespread interdisciplinary research, especially in mathematical modeling, to guide timely and practical intervention strategies. Ebola virus disease (EVD) is particularly challenging due to its high fatality rate, rapid symptom development, and potential for international spread. These factors make it not only a medical emergency but also a complex modeling problem that requires input from biology, sociology, and public health policy. Over the years, scientists have refined models to include not just disease transmission but also the societal, ecological, and decision-making contexts in which these outbreaks unfold. One of the earlier efforts to understand the spread of Ebola was made by Lekone and Finkenstädt.⁷ They used a stochastic SEIR model, which stands for Susceptible, Exposed, Infectious, and Recovered, to simulate how interventions like isolating patients and tracing their contacts could slow down the spread. What made their approach notable was how they factored in the randomness of real-world disease spread, something that becomes even more important in smaller populations where each case can dramatically change the outbreak's trajectory. Soon after, Ndanguza and colleagues⁸ focused on data from the 1995 outbreak in the Democratic Republic of Congo. They used statistical tools to dig into what actually happened during that outbreak. This kind of retrospective analysis helped build more accurate models for future use by highlighting how local healthcare quality and contact patterns influenced disease transmission over time. In 2014, Nigeria's experience with Ebola became a model for successful containment. Fasina et al.⁹ built a model to understand how early diagnosis, rapid response, and effective isolation measures helped stop the virus from spreading widely. Their work showed that local readiness, like good communication and coordination between agencies, can dramatically change the outcome of an epidemic, even in highly populated cities. Agusto et al.^10,11 added a cultural layer to traditional epidemiological models. They looked at how local customs, such as burial practices, affected transmission. For instance, traditional funerals often involve close physical contact with the deceased, which can unknowingly contribute to further infections. Their model demonstrated that overlooking these behaviors could lead to significant underestimation of outbreak severity. This idea was supported by Weitz and Dushoff,¹² who explored how the virus can still be contagious even after the host has died, making safe burial practices a critical public health intervention. On a broader scale, Ivorra and colleagues¹³ developed the Be-CoDiS model, short for Between-Countries Disease Spread. This model wasn't just about one community or one country, it was about how people, goods, and even the virus itself move across borders. Their model helped global health agencies better predict which countries might be at risk next and what steps could be taken to prepare. Berge et al.¹⁴ focused on making models that are both simple and useful, especially for countries with limited data and resources. They showed that even basic models could provide valuable guidance if they focused on the right variables, like how long someone remains infectious or how quickly they can be isolated. These small adjustments in the model had a big impact on the accuracy of predictions. Rachah and Torres¹⁵ approached the problem using control theory. Their models helped answer questions like: "How much should we spend on vaccinations? When should we focus more on public education?" Their findings suggested that adjusting strategies over time, rather than sticking to a fixed plan, often leads to better health outcomes and a more efficient use of resources. Abo and Smith¹⁶ explored what happens when vaccines are either available or unavailable. They built a model that accounted for changes in vaccine supply and how public acceptance might vary. Their work showed that even short interruptions in vaccination campaigns can drastically increase risk, highlighting the importance of consistent public health messaging and resource allocation. At the community level, Marais et al.¹⁷ stressed the need for cultural sensitivity in designing infection control programs. Their research showed that people are more likely to follow health guidelines when local leaders and traditions are involved in planning and communication. For example, religious leaders could help shape messages in a way that respects local customs while still promoting safe behavior. On the ecological side, Leroy et al.¹⁸ discovered that fruit bats likely serve as a natural reservoir for the Ebola virus. This means the virus can survive in these animals without making them sick, only to spill over to humans under certain conditions. Their work reinforced the idea that human encroachment into wildlife habitats, like forests, isn't just an environmental issue; it's a public health risk too. Burkett-Cadena and Vittor¹⁹ expanded on this by modeling how environmental changes, like deforestation, bring humans and wildlife into closer contact. Their models demonstrated how clearing forests or altering ecosystems can create new opportunities for viruses to jump from animals to people, making future outbreaks more likely. In terms of data-driven modeling, Kucharski et al.²⁰ used real-time data from Sierra Leone to adjust their models as the epidemic unfolded. This approach helped public health officials make quicker and better-informed decisions, such as when to deploy more resources or change intervention strategies. Their method marked a shift toward more responsive, real-time modeling that adapts to the outbreak as it evolves. Nishiura and Chowell²¹ brought all these ideas together in a comprehensive review. They evaluated dozens of different models to see what worked and what didn't. One of their key takeaways was that modeling accuracy improves dramatically when uncertainty, cultural diversity, and differences in healthcare access are taken seriously. They also called for better data-sharing practices to help researchers and governments act more cohesively during crises. Altogether, these studies reflect just how multifaceted Ebola modeling has become. Today's models aim to do more than just predict numbers, they try to capture the human, social, and ecological realities behind those numbers. Whether it's the behavior of a virus, a person, or a government, integrating all these layers gives us a better shot at preventing the next outbreak. The evolution of these modeling approaches shows that science and society must work hand-in-hand to effectively combat emerging infectious diseases. Agbomola and Loyinmi^22,23 analyzed the impact of human-bat interactions and control strategies on Ebola transmission using optimal control models, emphasizing the role of vector management and safe burial practices. Legrand et al.²⁴ provided critical insights into the effectiveness of public health interventions such as isolation, contact tracing, and the promotion of safe burial practices, which are crucial in breaking chains of transmission during outbreaks. Pandey et al.²⁵ explored the combined effects of vaccination campaigns, public awareness, and behavioral changes in reducing the spread of Ebola, offering a framework for understanding how early interventions could mitigate future outbreaks. Rivers et al.^23,26 demonstrated the power of real-time modeling in predicting outbreak trajectories, providing policymakers with tools to allocate resources effectively during the West African Ebola epidemic. Similarly, Althaus^27,28 employed dynamic transmission models to estimate key epidemiological parameters, such as the basic reproduction number R₀, for different stages of the epidemic. These models offered vital insights into the rate of spread and the effectiveness of interventions over time. Chretien et al.²⁹ underscored the role of environmental and ecological factors in Ebola outbreaks, highlighting how climatic conditions and land-use patterns influence disease emergence and persistence. While traditional models like the SIR and its derivatives provided a robust framework, they initially lacked the capacity to account for delays inherent in biological and societal processes. Delays, such as the incubation period of a disease, the time lag in symptom onset, or delays in implementing public health interventions, are crucial to accurately capturing disease dynamics. Addressing this limitation, delay differential equations (DDEs) emerged as a powerful mathematical tool to incorporate time-dependent processes into epidemiological models. Early contributions by Cooke and van den Driessche³⁰ demonstrated the utility of integrating delays into compartmental models, allowing for more realistic representations of disease progression and intervention outcomes. These models could, for example, simulate the effects of latent periods or the time-dependent impact of treatment protocols. Epidemiological modeling has long relied on ordinary differential equations (ODEs) to represent the dynamic spread of infectious diseases. These models, typified by the classical SIR framework, offer mathematical simplicity and analytical tractability. However, they inherently assume that all changes in a population's disease status occur instantaneously, overlooking the biological reality that many processes in disease transmission involve time delays. For instance, incubation periods, delays in symptom onset, and the lag between infection and the implementation of control measures are all time-dependent factors that shape disease dynamics in fundamental ways.⁵ By treating such processes as instantaneous, ODE-based models risk oversimplifying complex epidemiological phenomena. Delay differential equations (DDEs) address this critical limitation by explicitly incorporating time lags into the modeling framework.^30–32 This refinement allows for the representation of delayed infection responses, delayed treatment outcomes, and other temporally structured processes, offering a more biologically realistic depiction of disease progression.

Beyond epidemiology, the mathematical foundation of DDEs is well-established in dynamical systems theory, with results on stability analysis, bifurcation theory, and solution behavior extensively developed by researchers such as Hale and Lunel³³ and Beretta and Kuang.³⁴ These advances provide rigorous tools for understanding threshold phenomena, long-term behavior, and the conditions under which diseases can be controlled or eradicated. The flexibility of DDEs also extends to computational modeling. Shampine and Thompson³⁵ have shown how modern numerical methods can efficiently solve DDE systems, enabling the simulation of complex scenarios with time delays, such as delayed response to interventions or the accumulation of immunity in a population over time.

In addition, while ODE models offer foundational insights, DDEs provide a more accurate, realistic, and mathematically robust framework for modeling infectious diseases.^36,37 The explicit incorporation of delays bridges the gap between theoretical constructs and the biological realities of disease transmission, offering policymakers and public health experts more reliable tools for predicting outbreaks and designing effective interventions.^30–32,38

2. Model formulation

We examine the SEIRR model presented by Ref. 39, whose structure is illustrated in the schematic diagram below Figure 1 and Parameters Description are contained in Table 1. The model focuses on the host compartment of the epidemiological system, where the total population $N$ is divided into distinct epidemiological classes:

• Susceptible ( $S$ ): Individuals at risk of infection
• Exposed ( $E$ ): Infected individuals not yet infectious
• Infectious ( $I$ ): Individuals capable of transmitting the disease

The model further distinguishes between two types of Recovered Individuals: Even after recovering from Ebola with medical treatment, patients may still carry traces of the virus in certain parts of their bodies. Research shows the virus can linger in places like semen, eye fluid, breast milk, and spinal fluid - sometimes for weeks or months after symptoms disappear. The subdivision of the Recovered Individuals is:

• Recovered in Isolation ( $R_{1}$ ): Individuals who have recovered but remain under medical isolation due to potential viral persistence in body fluids
• Recovered without Isolation ( $R_{2}$ ): Individuals who have recovered and are not under isolation protocols

2.1 Model without delay differential equation

(1)

\begin{matrix} \frac{dS}{dt} = Λ - μS - βSI \\ \frac{dE}{dt} = βSI - (μ + θ) E \\ \frac{dI}{dt} = θE + ρ R_{2} - (μ + δ + cτ + (1 - c) τ) I \\ \frac{d R_{1}}{dt} = cτI - μ R_{1} \\ \frac{d R_{2}}{dt} = (1 - c) τI - (μ + ρ) R_{2} \end{matrix}

With Initial condition

S (t) \geq 0, E (t) \geq 0, I (t) \geq 0, R_{1} (t) \geq 0, R_{2} (t) \geq 0

Where

S (t) + E (t) + I (t) + R_{1} (t) + R_{2} (t) = N

2.2 Positivity of the model

1. Theorem

All the solutions of the system (1) are non-negative for $t \geq 0$ .

Proof

Let

\hat{y} = sup {t > 0 : S (t) \geq 0, E (t) \geq 0, I (t) \geq 0, R (t) \geq 0, R_{2} (t) \geq 0, P_{E} (t) \geq θ}

then

\hat{y} > 0

From the first equation in (1)

\frac{dS}{dt} = Λ - μS - βSI

\frac{dS}{dt} + (μ + βI) S = Λ

IF = e^{\int (μ + βI) dt} = e^{\int (μ + βI) dt}

e^{\int (μ + βI) dt} \frac{dS}{dt} + e^{\int (μ + βI) dt} (μ + βI) S = Λ e^{\int (μ + βI) dt}

\frac{d}{dt} {S e^{\int (μ + βI) dt}} = Λ e^{\int (μ + βI) dt}

Let

(S (0) = S_{0}, E (0) = E_{0}, I (0) = I_{0}, R_{1} (0) = R_{1_{0}}, R_{2} (0) = R_{2_{0}})

S (t) = e^{- \int_{0}^{t} (μ + βI) dt} (S (0) + \int_{0}^{t} Λ e^{\int_{0}^{s} (μ + βI) dt} dt)

As $S (t) \geq 0$ , $S$ is non-negative.

Similarly, we prove that other compartments

E (t), I (t), R_{1} (t), R_{2} (t)

are positive

\forall t \geq 0 .

2.3 Boundedness of the model

2. Theorem

The analytic solutions of the system of equations in (1) are bounded.

Proof

By summing the system of equations in (1)

\frac{dS}{dt} + \frac{dE}{dt} + \frac{dI}{dt} + \frac{d R_{1}}{dt} + \frac{d R_{2}}{dt}

which simplifies to:

(2)

\frac{d (S + E + I + R_{1} + R_{2})}{dt} = Λ - μ (S + E + I + R_{1} + R_{2}) - δI

which gives:

\frac{d (S + E + I + R_{1} + R_{2})}{dt} \leq Λ - μ (S + E + I + R_{1} + R_{2})

Recall $S + E + I + R_{1} + R_{2} = N$

\frac{dN}{dt} \leq Λ - μN

Integrating:

N (t) \leq \frac{Λ}{μ} + C e^{- μt}

As $t \to \infty$ , $C e^{- μt}$ decays to zero:

N (t) \leq \frac{Λ}{μ}

Hence:

\underset{t \to \infty}{limsup} (S + E + I + R_{1} + R_{2}) \leq \frac{Λ}{μ}

The feasible region for the system of equations in (1) is:

Ω = {(S, E, I, R_{1}, R_{2}); S + E + I + R_{1} + R_{2} \leq \frac{Λ}{μ}, S \geq 0, E \geq 0, I \geq 0, R_{1} \geq 0, R_{2} \geq 0}

3. Basic Reproductive Number ( $R_{0}$ )

From the system of equation in (1)

(3)

\begin{matrix} \frac{dS}{dt} = Λ - μS - βSI \\ \frac{dE}{dt} = βSI - (μ + θ) E \\ \frac{dI}{dt} = θE + ρ R_{2} - (μ + δ + cτ + (1 - c) τ) I \\ \frac{d R_{1}}{dt} = cτI - μ R_{1} \\ \frac{d R_{2}}{dt} = (1 - c) τI - (μ + ρ) R_{2} \end{matrix}

At DFE

[S, E, I, R_{1}, R_{2}] = [\frac{Λ}{μ}, 0, 0, 0, 0]

Let X^{'} = {[E, I, R_{2}]}^{T} (Disease states of infection)

X^{*} = F (X^{*}) - V (X^{*})

F = [\begin{matrix} 0 & βSI & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}]

V = [\begin{matrix} (μ + θ) E & 0 & 0 \\ - θE & (μ + δ + c_{1} τ + (1 - c) τ) I & - ρ R_{2} \\ 0 & - (1 - c) τI & (μ + ρ) R_{2} \end{matrix}]

{F (x) = \frac{\partial F}{\partial x} |}_{at DFE}

{V (x) = \frac{\partial V}{\partial x} |}_{at DFE}

F (x) = [\begin{matrix} 0 & βS & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}] = [\begin{matrix} 0 & β \frac{Λ}{μ} & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix}]

V (x) = [\begin{matrix} (μ + θ) & 0 & 0 \\ 0 & (μ + δ + c_{1} τ + (1 - c) τ) & - ρ \\ 0 & - (1 - c) τ & (μ + p) \end{matrix}]

{FV}^{- 1} = (\begin{matrix} \frac{β Λ θ (μ + ρ)}{μ (θ + μ) ((μ + ρ) (c τ + δ + μ) - (c - 1) μτ)} & \frac{β Λ (μ + ρ)}{μ ((μ + ρ) (c τ + δ + μ) - (c - 1) μτ)} & \frac{β Λ ρ}{μ ((μ + ρ) (c τ + δ + μ) - (c - 1) μτ)} \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{matrix})

R_{0} = \frac{β Λ θ (μ + ρ)}{μ (θ + μ) ((μ + ρ) (cτ + δ + μ) - (c - 1) μτ)}

4. Model with delay differential equation

It is common to observe that infected individuals do not transmit the infection to others immediately. There is a time delay before any confirmed case can occur after the exposure. This delay is known as the latent period. This introduces a gap between exposure and when the infection becomes transmissible. In mathematical form, this latent period is modeled as a time delay in the infection process. Time delays can destabilize the equilibrium point of the system, leading to a Hopf bifurcation. This results in periodic oscillations. This motivates the inclusion of the delay in the model reflecting that individuals infected at time $t - τ_{2}$ can spend at time $t$ .

Let’s take into consideration a delayed SIRR model with a saturation incidence rate and exponential birth rate:

\frac{S (t - τ_{2}) I (t - τ_{2}) e^{- μ_{2} τ_{2}}}{1 + mI (t - τ_{2})}

where $m$ is a constant related to the saturation effect.

The system of equations below represents the delay differential equations from this model:

(4)

\begin{matrix} \frac{dS}{dt} = Λ - μS - \frac{βS (t - τ_{2}) I (t - τ_{2}) e^{- μ_{2} τ_{2}}}{1 + mI (t - τ_{2})} \\ \frac{dI}{dt} = \frac{βS (t - τ_{2}) I (t - τ_{2}) e^{- μ_{2} τ_{2}}}{1 + mI (t - τ_{2})} + ρ R_{2} - (μ + δ + cτ + (1 - c) τ) I \\ \frac{d R_{1}}{dt} = cτI - μ R_{1} \\ \frac{d R_{2}}{dt} = (1 - c) τI - (μ + ρ) R_{2} \end{matrix}

From the system of equations, we have in (4)

We can see that $R_{1}$ and $R_{2}$ compartments are independent of $S$ and $I$ .

Thus, it’s enough to consider the following reduced system of equations for our study:

(5)

\begin{matrix} \frac{dS}{dt} = Λ - μS - \frac{βS (t - τ_{2}) I (t - τ_{2}) e^{- μ_{2} τ_{2}}}{1 + mI (t - τ_{2})} \\ \frac{dI}{dt} = \frac{βS (t - τ_{2}) I (t - τ_{2}) e^{- μ τ_{2}}}{1 + mI (t - τ_{2})} - (μ + δ + cτ + (1 - c) τ) I \end{matrix}

with initial conditions

S (t) \geq 0, I (t) \geq 0

.

4.1 Equilibrium and stability analysis

The equilibria of system (5) are found by setting the right-hand side of the system to zero. A stable endemic equilibrium point is a fixed point in a system where if the system is slightly disturbed, it will eventually return to that point. Therefore, the equilibrium solutions of a system with a time delay are the same as those of the corresponding system without delay.

4.2 Disease Free Equilibrium

The D.F.E can be found by setting the system of equation in $3$ to zero. Which is $P = {\frac{Λ}{μ}, 0}$ i.e. $P = {S_{0}, 0}$

The basic reproductive number $R_{0}$ : This is the number of secondary infections caused by a single infected individual in a completely susceptible population.

From system (4), we compute $R_{0}$ by setting

\frac{dI}{dt} \geq 0, assuming there is no delay I (t) \to 0 .

(6)

\frac{dI}{dt} = βSI + ρ R_{2} - (μ + δ + cτ + (1 - c) τ) I \geq 0

Assuming $R_{2} \to 0$ , we have:

β S_{0} \geq (μ + δ + cτ + (1 - c) τ) \Rightarrow \frac{β S_{0}}{μ + δ + cτ + (1 - c) τ} \geq 1

Thus,

\frac{β Λ}{μ (μ + δ + cτ + (1 - c) τ)} \geq 1

Finally, we get:

R_{0} = \frac{β Λ}{μ (μ + δ + cτ + (1 - c) τ)}

3. Theorem

The disease-free equilibrium $P$ is locally asymptotically stable if $R_{0} < 1$ and unstable if $R_{0} > 1$ $\forall τ_{2} \geq 0$

Proof

Let $F$ and $V$ be $3 \times 3$ matrices representing the infection and transition dynamics, respectively, evaluated at the DFE. Assume that $F \geq 0$ (non-negative matrix) and $V$ is a nonsingular $M$ -matrix (i.e., a matrix whose inverse is non-negative and has eigenvalues with positive real parts).

J = F - V = (\begin{matrix} F_{11} - V_{11} & F_{12} - V_{12} & F_{13} - V_{13} \\ F_{21} - V_{21} & F_{22} - V_{22} & F_{23} - V_{23} \\ F_{31} - V_{31} & F_{32} - V_{32} & F_{33} - V_{33} \end{matrix})

The basic reproductive number

R_{0} = ρ (F V^{- 1})

CASE 1

When R $_{0}$ $<$ 1

We know that $F \geq 0$ and $V$ is a non-singular $M$ -matrix. Hence, $R_{0} = ρ (F V^{- 1}) > 1$ , and all eigenvalues of $(F - V)$ have negative real parts.

Since $R_{0} = ρ (F V^{- 1}) > 1$ , we have that:

(I - F V^{- 1}) \geq 0

If $(V - F) \geq V (I - V^{- 1} F)$ , then ${(V - F)}^{- 1} \geq {(I - V^{- 1} F)}^{- 1}$ $V^{- 1}$ . Recall that $V^{- 1} \geq 0$ .

It follows that:

{(V - F)}^{- 1} \geq 0

Therefore, $V - F$ is a non-singular $M$ -matrix, and by the characteristics of non-singular $M$ -matrices, all the eigenvalues of $V - F$ have positive real parts.

Consequently, all eigenvalues of $F - V$ have negative real parts, implying that the DFE is locally asymptotically stable.

CASE 2

When $R_{0} > 1$

We can prove that if $R_{0} > 1$ , the DFE is unstable. Assume $R_{0} > 1$ and suppose by contradiction that all eigenvalues of $F - V$ have non-negative real parts.

λ : (F - V) \geq 0, \forall i, where i \in ℕ

Since $ϵ > 0$ is arbitrary:

(I + ϵI) - F V^{- 1}

For any $ϵ > 0$ , this matrix $(I + ϵI) - F V^{- 1}$ remains a non-singular $M$ -matrix. Hence, it has all positive eigenvalues.

Thus, when $R_{0} > 1$ , the matrix $F - V$ has at least one eigenvalue with a positive real part, which implies the DFE is unstable.

CASE 3

When $R_{0} = 1$

When $R_{0} = 1$ we can see that $β = β^{*} = \frac{μ (μ + δ + cτ + (1 - c) τ)}{Λ}$

When $β^{*}$ is the bifurcation parameter, hence the system has a zero eigenvalue (simple eigenvalue) and the other is a negative eigenvalue, i.e. nonhyperbolic equilibrium.

To study the thermal stability of the nonhyperbolic equilibrium, we will need to apply the center manifold theory [Sastry(1999)]⁴² and Theorem 4.1 of [castillo-chaves and song (2004)]⁴⁰

To describe this, let’s consider the system a general system for ordinary differential equations.

\frac{dx}{dt} = f (x, β^{*}), f : ℝ^{n} \times ℝ \to ℝ^{n} and f \in C^{2} (ℝ^{n} \times ℝ)

Let $x_{1} = S, x_{2} = I$ , the system in (5) becomes:

(7)

\begin{matrix} F_{1} = \frac{d x_{1}}{dt} = Λ - μ x_{1} - \frac{β^{*} x_{1} (t - τ) x_{2} (t - τ) e^{- μ_{2} τ_{2}}}{1 + m x_{2} (t - τ)} \\ F_{2} = \frac{d x_{2}}{dt} = \frac{β^{*} x_{1} (t - τ) x_{2} (t - τ) e^{- μ_{2} τ_{2}}}{1 + m x_{2} (t - τ)} + ρ x_{4} - (μ + δ + cτ + (1 - c) τ) x_{2} \end{matrix}

At $R_{0} = 1, β = β^{*} = \frac{μ (μ + δ + cτ + (1 - c) τ)}{Λ}$

at $DFE (x_{1} = \frac{Λ}{μ}, x_{2} = 0)$

The jacobian matrix of the system above around $DFE$

Det = [\begin{matrix} - μ & \frac{- β^{*} Λ}{μ} \\ 0 & 0 \end{matrix}]

Let’s consider the right eigenvector associated with 0 eigenvalue which is $u = [u_{1}, u_{2}] = [\frac{- β Λ}{μ^{2}}, 1]$ and left eigen vector associated with the zero eigenvalue (simple eigenvalue is $w = [w_{1}, w_{2}] = [0, 1]$ )

When we consider the partial derivative associated with the system.

\begin{matrix} \frac{\partial^{2} F_{1}}{\partial x_{1} \partial x_{2}} = - β^{*}, \frac{\partial^{2} F_{1}}{\partial x_{1}^{2}} = 0 \\ \frac{\partial^{2} F_{2}}{\partial x_{1} \partial x_{2}} = β^{*}, \frac{\partial^{2} F_{2}}{\partial x_{2}^{2}} = \frac{- 2 Λ m β^{*}}{μ} \\ \frac{\partial^{2} f_{2}}{\partial x_{2} \partial β^{*}} = \frac{Λ}{μ} \end{matrix}

Now we can apply Theorem 4.1[Castillo-chaves and song[2004]⁴⁰

When the formula for $a$ and $b$ are

\begin{matrix} a = \sum_{k = 1, j = 1}^{n} w_{k} u_{i} u_{j} \frac{\partial^{2} f_{k}}{\partial x_{i} \partial x_{j}} (0, 0) \\ b = \sum_{k = 1, j = 1}^{n} w_{k} u_{i} \frac{\partial^{2} k}{\partial x_{1} \partial β} (0, 0) \\ a = \sum_{k = 1}^{n} \sum_{1 = 1}^{n} \sum_{j = 1} w_{k} u_{i} u_{j} \frac{\partial^{2} F_{2}}{\partial x \partial y} (0, 0) . \\ = - (\frac{2 β^{*} Λ (β^{*} + μm)}{μ^{2}}) < 0 \end{matrix}

b = u_{2} (\frac{Λ}{μ}) = \frac{Λ}{μ} > 0

1. Remarks:

From Theorem 4.1[Castillo chaves and song [2004].⁴⁰ Since $a < 0$ and $b > 0,$ the system undergoes a transcritical bifurcation at $R_{0} = 1$ as the bifurcation parameter is varied.

4.3 Endemic equilibrium

To establish the conditions for the existence of endemic equilibrium $P^{*} = (S^{*}, I^{*})$ , when the disease remains present in the population, the system of equation (5) is reformulated to derive the fixed points for $S^{*}$ and $I^{*}$ yielding the following results;

(1) $I^{*} = \frac{Λ - μ s^{*}}{μ s^{*} m + β s^{*} - Λ m}$
(2) ${s^{*}}$

$S^{*}$ is given by the quadratic equation;

h_{1} s^{2 *} + h_{2} s^{*} + h_{3} = 0

Where

\begin{matrix} h_{1} = zμ - μβ - μ^{2} m \\ h_{2} = Λ μm + Λ β + μm + ρ R_{2} Λ m + ρ R_{2} β \\ h_{3} = - Λ^{2} m - ρ R_{2} Λ m - z Λ \end{matrix}

Where $z = (μ + δ + cτ + (1 - c) τ)$

By Rene Descarte’s rule of sign we can say that the quadratic equation have a unique positive real root if the following condition holds:

h_{1} > 0, h_{2} > 0, h_{3} < 0

Local stability of the endemic equilibrium $P^{*}$ is analysed as follows; The characteristic equation of system evaluated at $P^{*}$

The characteristic equation

λ^{2} + J_{1} λ + k_{1} + (J_{2} λ + k_{2}) e^{- λ τ_{2}}

When we have

\begin{matrix} J_{1} = μ + z + \frac{β I^{*}}{(1 + m I^{*})}, k_{2} = \frac{mβ s^{*}}{{(1 + m I^{*})}^{2}} \\ k_{1} = zμ + \frac{zβ I^{*}}{(1 + m I^{*})} \\ J_{2} = \frac{- β s^{*}}{{(1 + m I^{*})}^{2}} \end{matrix}

To ensure that the endemic equilibrium $P^{*}$ is locally asymptotically stable, we need both eigenvalues $λ$ to have negative real parts.

According to the routh-Hurwitz criterion, this happens if and only if;

\begin{matrix} J > 0 \\ k > 0 \end{matrix}

Hence the theorem:

4. Theorem

for $τ_{2} = 0$ the endemic equilibrium $P^{*} = (S^{*}, I^{*})$ is locally asymptotically stable if both of the following conditions hold simultaneously;

(1) $\frac{s^{*}}{I^{*}} \leq \frac{μ}{z}$
(2) $\frac{s^{*}}{I^{*}} \leq I$
$\begin{matrix} J = J_{1} + J_{2} = μ + z + \frac{β I^{*}}{(1 + mI)} - \frac{β s^{*}}{{(1 + mI)}^{2}} \\ = μ + z + \frac{βI + βm I^{2 *} - β s^{*}}{{(1 + mI)}^{2}} > 0 \\ k = k_{1} + k_{2} = zμ + \frac{zβ I^{*}}{(1 + m I^{*})} + \frac{μβ s^{*}}{{(1 + m I^{*})}^{2}} > 0 \end{matrix}$

Hence $J, k > 0$ and satisfies Routh-Hurwitz criterion.

5. Hopf Bifurcation analysis

The characteristic equation is given by:

(8)

λ^{2} + J_{1} λ + k_{1} + (J_{2} λ + k_{2}) e^{- λ τ_{2}} = 0 .

Substituting $λ = iw$ :

(9)

(- w^{2} + iw J_{1} + k_{1}) + (i J_{2} + k_{2}) (cos (w τ_{2}) - i sin (w τ_{2})) = 0 .

Separating real and imaginary parts:

\begin{matrix} Real part : - w^{2} + k_{1} + k_{2} cos (w τ_{2}) - w J_{2} sin (w τ_{2}) = 0, \\ Imaginary part : w J_{1} + w J_{2} cos (w τ_{2}) - k_{2} sin (w τ_{2}) = 0 . \end{matrix}

Squaring and adding these parts eliminates $cos (w τ_{2})$ and $sin (w τ_{2})$ , yielding:

(10)

w^{4} + (J_{1}^{2} + J_{2}^{2} - 2 k_{1}) w^{2} + (k_{1}^{2} - k_{2}^{2}) = 0 .

Let $ξ_{1} = w^{2}$ . The equation becomes:

(11)

ξ_{1}^{2} + P ξ_{1} + Δ = 0,

where:

P = J_{1}^{2} + J_{2}^{2} - 2 k_{1}, Δ = k_{1}^{2} - k_{2}^{2} .

then we get,

ξ_{1}^{2} + P ξ_{1} + Δ = 0 .

From the equation above

τ_{n} = \frac{1}{w} arccos (\frac{(k_{2} - J_{1} J_{2}) w^{2} - k_{1} k_{2}}{J_{2} w^{2} + k_{1}^{2}}) + \frac{2 nπ}{w}, n = 0, 1, 2, \dots

Assuming the discriminant $P^{2} - 4 Δ < 0$ .

λ^{2} + J_{1} λ + k_{1} + (J_{2} λ + k_{2}) e^{- λ τ_{2}} = 0

Differentiate with respect to $τ_{2}$

(12)

2 λ \frac{dλ}{d τ_{2}} + J_{1} \frac{dλ}{d τ_{2}} + (J_{2} \frac{dλ}{d τ_{2}} e^{- λ τ_{2}} + (J_{2} λ + k_{2}) \frac{d}{d τ_{2}} (e^{- λ τ_{2}})) = 0

\frac{d}{d τ_{2}} (e^{- λ τ_{2}}) = - λ e^{- λ τ_{2}} - τ_{2} \frac{dλ}{d τ_{2}} e^{- λ τ_{2}}

2 λ \frac{dλ}{d τ_{2}} + J_{1} \frac{dλ}{d τ_{2}} + J_{2} \frac{dλ}{d τ_{2}} e^{- λ τ_{2}} + (J_{2} λ + k_{2}) (- λ e^{- λ τ_{2}} - τ_{2} \frac{dλ}{d τ_{2}} e^{- λ τ_{2}}) = 0

2 λ \frac{dλ}{d τ_{2}} + J_{1} \frac{dλ}{d τ_{2}} + J_{2} \frac{dλ}{d τ_{2}} e^{- λ τ_{2}} - (J_{2} λ + k_{2}) λ e^{- λ τ_{2}} - (J_{2} λ + k_{2}) τ_{2} \frac{dλ}{d τ_{2}} e^{- λ τ_{2}} = 0

(2 λ + J_{1} + J_{2} e^{- λ τ_{2}} - (J_{2} λ + k_{2}) τ_{2} e^{- λ τ_{2}}) \frac{dλ}{d τ_{2}} = (J_{2} λ + k_{2}) λ e^{- λ τ_{2}}

\frac{dλ}{d τ_{2}} = \frac{(J_{2} λ + k_{2}) λ e^{- λ τ_{2}}}{2 λ + J_{1} + J_{2} e^{- λ τ_{2}} - (J_{2} λ + k_{2}) τ_{2} e^{- λ τ_{2}}}

{[\frac{dλ}{d τ_{2}}]}^{- 1} = \frac{d τ_{2}}{dλ} = \frac{2 λ + J_{1} + J_{2} e^{- λ τ_{2}} - (J_{2} λ + k_{2}) τ_{2} e^{- λ τ_{2}}}{(J_{2} λ + k_{2}) λ e^{- λ τ_{2}}}

Recall that:

\begin{matrix} λ^{2} + J_{1} λ + k_{1} + (J_{2} λ + k_{2}) e^{- λ τ_{2}} \\ - λ^{2} + J_{1} λ + k_{1} = (J_{2} λ + k_{2}) e^{- λ τ_{2}} \end{matrix}

Replacing this in the equation.

\begin{matrix} \frac{2 λ + J_{1}}{- λ (λ^{2} + J_{1} λ + k_{1})} + \frac{J_{2}}{λ (J_{2} λ + k_{2})} - \frac{τ_{2}}{λ} \\ {\frac{d}{d τ_{2}} {(Re (λ))}^{- 1} = Re {(\frac{dλ}{d τ_{2}})}^{- 1} |}_{λ = iρ} \\ \begin{matrix} = \frac{1}{w} (\frac{2 w (w^{2} - k_{0}) + J_{1} w}{{(w^{2} - k_{0})}^{2} + {(J_{1} w)}^{2}} - \frac{J_{1} w}{{(J_{1} w)}^{2} + k_{2}^{2}}) \\ (w^{2} - k_{1}) + {(J_{1} w^{2})}^{2} = {(J_{2} w)}^{2} + k_{2}^{2} \end{matrix} \end{matrix}

Given the condition $J_{1}^{2} + J_{2}^{2} - 2 k_{1}$ $> 0$ , it follows that:

{\frac{d}{dτ} Re {(λ)}^{- 1} |}_{λ = iw} > 0,

which ensures that the transversality condition is satisfied. This, in turn, leads to the occurrence of a Hopf bifurcation at

ω = w

and

τ_{2} = τ (0)

. As a result, the equilibrium

P^{*}

of system (5) is asymptotically stable for

τ_{2} \in (0, τ (0)

, and undergoes a Hopf bifurcation at

τ_{2} = τ (0)

, marking a significant change in the system’s dynamical behavior.

6. Global asymptotically stability

From equation (5), we have

\frac{dS}{dt} = Λ - μS - \frac{βS (t - τ_{2}) I (t - τ_{2}) e^{- μ τ_{2}}}{1 + mI (t - τ_{2})}

\frac{dI}{dt} = \frac{βS (t - τ_{2}) I (t - τ_{2}) e^{- μ τ_{2}}}{1 + mI (t - τ_{2})} - ρ R_{1} + (μ + δ + cτ + (1 - c) τ) I

can be re-written as

\frac{dS}{dt} = Λ - μS - U (S (t), I (t))

\frac{dI}{dt} = U (S (t), I (t)) + V (I (t))

The Lyapunov direct method, as formulated by Huang et al,⁴¹ provides a systematic approach for establishing the global stability of equilibria in dynamical systems. This method involves the construction of a suitable scalar-valued function, known as a Lyapunov function, which satisfies specific positivity and derivative conditions along system trajectories. In this section, we construct an appropriate Lyapunov functional adapted to the structure of the model (5) and employ it to thoroughly analyze the global stability of the system’s equilibrium points.

The incidence function $U (S (t), I (t))$ is expressed as the product $U (S (t)) W (I (t))$ .

V (S (t)) = βS (t)

W (I (t)) = \frac{I (t)}{1 + mI (t)}

Y (I (t)) = ρ R_{1} + zI, where z = (μ + δ + cτ (1 - c) τ)

and impose the following assumptions:

1. The functions $U$ and $W$ satisfy $U (0) = W (0) = 0$ and their first derivatives $V^{'} (S)$ and $W^{'} (I)$ are strictly positive for all $S, I > 0$ .
2. Either $\frac{W (I (0))}{W (I^{*})} \geq 1$ or $\frac{Y (I (t))}{V (S (t)) W (I (t))} \leq 1$ for all $S, I \geq 0$ .
3. The function $W$ is concave in $I$ , i.e., $W^{'} (I) > 0$ and $\frac{\partial^{2} W (I)}{\partial I^{2}}$ for all $S, I > 0$ .

On this premise, we state the following theorem:

5. Theorem

1. Under the conditions 1 and 2, the equilibrium point $P^{*} (S^{*}, I^{*})$ is globally stable for any delay parameter $τ_{2} \geq 0$ .
2. If conditions (1)–(3) hold and $R_{0} < 1$ , the DFE is stable for $τ_{2} \geq 0$ .

Proof

If we define

F_{1} (S (t), I (t)) = S (t) - V (S (t)) \int_{S^{*}}^{S (t)} \frac{dϕ}{V (ϕ)} + I (t) - W (I (t)) \int_{I^{*}}^{I (t)} \frac{dϕ}{W (ϕ)}

Differentiating by Leibniz’s rule,

\frac{\partial F_{1}}{\partial S} = 1 - \frac{V (S^{*})}{V (S)}

\frac{\partial F_{1}}{\partial I} = 1 - \frac{W (I^{*})}{W (I)}

Since

\frac{\partial^{2} F_{1}}{\partial S^{2}} = \frac{V (S) U (S^{*})}{{(V (S))}^{2}} > 0 and

\frac{\partial^{2} F_{1}}{\partial I^{2}} = \frac{W (I) U (I^{*})}{{(W (I))}^{2}} > 0

\frac{\partial^{2} F_{1}}{\partial S \partial I} = 0

Thus, $P^{*} (S^{*}, I^{*})$ is globally minimum.

Taking the total derivatives

\frac{d F_{1}}{dt} = \frac{dS}{dt} [1 - \frac{V (S^{*})}{V (S)}] + \frac{dI}{dt} [1 - \frac{W (I^{*})}{W (I)}]

\frac{d F_{1}}{dt} = [V (S^{*}) W (I^{*}) - μ S^{*} - \frac{V (S (t - τ_{2}) W (I (t - τ_{2})))}{1 + mI (t - τ_{2})}]

\times [1 - \frac{V (S, x^{*})}{V (S)}] + [\frac{V (S (t - τ_{2})) W (t - τ_{2})}{1 + WI (t - 2 ϕ)} + ρ R_{1} + ZI] [1 - \frac{W (I^{*})}{W (I)}]

Let’s assume

For $F_{2}$ :

\begin{matrix} F_{2} = \int_{0}^{τ_{2}} \frac{W (I (t - ϕ))}{W (t^{*})} - 1 - ln [\frac{W (I (t - ϕ))}{W (t)}] dϕ \\ \frac{d F_{2}}{dt} = \frac{d}{dt} \int_{0}^{τ_{2}} \frac{W (I (t - ϕ))}{W (I (t^{*}))} - 1 - ln [\frac{W (I (t - ϕ))}{W (I (t))}] dϕ \\ = - \int_{0}^{τ_{2}} \frac{d}{dϕ} [\frac{W (I (t \int ϕ))}{W (I^{*})} \int 1 - ln (\frac{W (I (t - ϕ))}{W (I (t))})] dϕ \\ - \frac{W (t - τ_{2})}{W (I)} + \frac{W (I (t))}{W (I^{*})} + ln (\frac{W (I (t - τ_{2}))}{W (I (t))}) \end{matrix}

Let’s define the Lyapunov function

F = F_{1} + G F_{2} where G = V (S^{*}) W (I^{*})

\frac{dF}{dt} = \frac{d F_{1}^{*}}{dt} + G \frac{d F_{2}}{dt}

\begin{matrix} \frac{dF}{dt} = [V (S^{*}) W (I^{*}) - μ S^{*} - \frac{V (S (t - τ_{2})) W (I (t - τ_{2})))}{1 + mI (t - τ_{2})}] [1 - \frac{V (S^{*})}{V (S)}] \\ + [\frac{V (S ((t - τ_{2}) W (t - τ_{2})))}{1 + mI (t - τ_{2})} + y (I (t))] [\frac{1 - W (I (t))}{W (I (t))}] \\ + V (S^{*}) W (T^{*}) [\frac{- W (I (t - τ_{2}))}{W (I^{*})} + \frac{W (I (T))}{W (I^{*})} + ln \frac{WI (t - τ_{2})}{WI (t)}] \end{matrix}

Let’s use;

\begin{matrix} ln \frac{W (I (t - τ_{2}))}{W (I (t))} = ln (\frac{V (S^{*} (1 + ln (t - τ_{2})))}{V (S (t - τ_{2}))}) + ln [\frac{V (S (t - τ_{2}) W (I (t - τ_{2}))}{W (I (t) V (S^{*} (1 + ln (t - τ_{2})))}] \\ \begin{matrix} \frac{dF}{dt} = [1 - \frac{V (S *)}{V (S)} + ln \frac{V (S^{*}) (1 + mI (t - τ_{2}))}{V (s (t - τ_{2}))}] V (S^{*}) W (I^{*}) \\ + [V (S^{*}) - V (S))] [μ S^{*} + \frac{V (S (t - τ_{2})) W (I (t - τ_{2}))}{1 + mI (t - τ_{2})}] \frac{1}{V (S)} \\ \begin{matrix} + [V (S^{*}) W (I^{*})] [1 - \frac{V (S (t - τ_{2})) W (I (t - τ_{2})}{V (S^{*}) (1 + mI (t - τ_{2}))}] + ln [\frac{V (S (t - τ_{2})) W (I (t - τ_{2}))}{W (I (t) V (S^{*})) (1 + mI (t - τ_{2}))}] \\ + V (S^{*}) W (I^{*}) [1 - \frac{W (I (t))}{W (I^{*})}] [1 - \frac{W (I^{*})}{W (I)}] \end{matrix} \end{matrix} \end{matrix}

If we have

[1 - \frac{V (S^{*})}{V (S)} + ln (\frac{V (S^{*}) (1 + mI (t - τ_{2}))}{V (S (t - τ_{2})}))] < 0

and

[1 - \frac{V (S (t - τ_{2})) W (I (t - τ_{2}))}{V (S^{*}) (1 + mI (t - τ_{2}))} + μ \frac{V (S (t - τ_{2})) W (I (t - τ_{2}))}{W (I (t)) V (S^{*}) (1 + mI (t - τ_{2}))}] < 0

Also

V (S) > V (S^{*}) \forall S (t) > S^{*},

hence

[V (S^{*}) - V (S)] [μ S^{*} + \frac{V (S (t - τ_{2})) W (I (t - τ_{2}))}{1 + mI (t - τ_{2})}] \leq 0

Hence we can say that,

\frac{dF}{dt} \leq 0,

This is sufficient for Corollary 5.2 of [Kuang (1993).³¹ We can then say that the endemic state $P^{*} (S^{*}, I^{*})$ is globally asymptotically stable.

7. Discussion of Results

1. Figures (2)-(8) and Figure 17: Dynamics Comparison (With and Without Delay) [ Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 , Figure 8.]
The dynamics comparison plot demonstrates the effect of incorporating a delay ( $τ_{2}$ ) in the recovery process on the population dynamics of susceptible ( $S$ ) and infected ( $I$ ) individuals.
- • Without Delay:
  The system without delay ( $τ_{2} = 0$ ) exhibits a faster approach to equilibrium, reflecting a more straightforward trajectory to either disease-free equilibrium (DFE) or endemic equilibrium (EE), depending on parameter values. This is consistent with classical models, where recovery and transitions between compartments occur instantaneously or without delays.
- • With Delay:
  Introducing a delay ( $τ_{2} > 0$ ) introduces oscillatory behavior in the infected population due to the feedback loop created by the delayed transitions. These oscillations may represent cycles of outbreaks and remissions, commonly observed in zoonotic diseases like Ebola, where external factors such as virus persistence in bodily fluids or immune system delays affect recovery rates.
  The presence of delay slows down convergence to equilibrium and may indicate higher disease prevalence over time. These dynamics underscore the critical role of delayed feedback in influencing outbreak persistence and stability in real-world scenarios.
2. Figure 12 and Figures 13–Figure 15: Phase Space (Trajectories in $S, I, R_{1}$ and $S, I, R_{2}$ )
The phase space plot provides insights into the long-term behavior of the model with and without delay.
- • Disease-Free Equilibrium (DFE):
  The DFE point (marked as black in the plot) corresponds to the scenario where the infection dies out completely ( $I = 0$ , $R_{1} = R_{2} = 0$ ). This equilibrium is stable for small transmission rates ( $β$ ) or when the basic reproduction number ( $R_{0}$ ) is less than 1.
- • Endemic Equilibrium (EE):
  The endemic equilibrium (marked as green) represents a state where the infection persists in the population. This occurs when $R_{0} > 1$ , allowing for a sustained infected population. For the delayed system, trajectories often show oscillations before converging to the EE, reflecting periodic outbreaks.
  The trajectory with delay exhibits more complex behavior (such as spiraling orbits) before stabilization, highlighting the destabilizing impact of delayed recovery on the population dynamics.
3. Figure 16: Bifurcation Analysis (Peak Infection vs. Delay $τ_{2}$ )
The bifurcation diagram investigates how varying the delay ( $τ_{2}$ ) influences peak infection levels.
- • Without Delay ( $τ_{2} = 0$ ):
  Peak infection levels remain relatively constant as the system lacks the feedback loop induced by delay. The dynamics stabilize quickly at the endemic equilibrium.
- • With Delay ( $τ_{2} > 0$ ):
  Introducing delay leads to increasing peak infection levels as $τ_{2}$ grows. This is a direct consequence of the prolonged infectious period caused by delays in transitions to recovery ( $R_{1}$ and $R_{2}$ ). Moreover, as $τ_{2}$ increases beyond a critical threshold, the system exhibits oscillatory outbreaks, as seen in the dynamics plot.
  The bifurcation plot effectively distinguishes between the disease-free region ( $R_{0} < 1$ ) and the endemic region ( $R_{0} > 1$ ) by showing when peak infection levels transition from zero to non-zero values. The critical value of $τ_{2}$ at which this transition occurs aligns with thresholds derived from the model’s reproduction number.
4. Isolation Parameter $c$
The isolation parameter $c$ reflects the fraction of recovered individuals ( $R_{1}$ ) who effectively isolate themselves and avoid further transmission. The two plots ( Figure 9 and Figure 10) show how $R_{0}$ depends on $c$ , and the behavior is analyzed for $R_{0} > 1$ and $R_{0} < 1$ :
- • Figure 9 ( $R_{0} > 1$ ):
  With the range $[4.3, 5.9]$ , there is a decreasing relationship between $c$ and $R_{0}$ , indicating that as the isolation parameter $c$ increases, $R_{0}$ decreases. Increasing $c$ reduces the effective transmission by isolating recovered individuals, helping to lower $R_{0}$ .
- • Figure 10. ( $R_{0} < 1$ ):
  With the range $[0.43, 0.59]$ , there is a decreasing relationship between $c$ and $R_{0}$ , but $R_{0}$ values are already below 1. Increasing $c$ further accelerates the decline, making it harder for the infection to persist in the population. Here, lower isolation efforts (small $c$ ) are sufficient to maintain $R_{0} < 1$ .
5. Figure 11: Sensitivity Analysis of $R_{0}$
Figure 11 presents the sensitivity analysis of the basic reproduction number ( $R_{0}$ ) with respect to model parameters. The results reveal that the transmission rate ( $β$ ), recruitment rate ( $Λ$ ), and recovery rate into $R_{2}$ ( $ρ$ ) are the parameters that contribute most significantly to $R_{0}$ .

Figure 1. Schematic diagram.

Table 1. Parameters of the SEIRR Ebola transmission model.

Parameter	Description
$Λ$	Recruitment rate of susceptible individuals
$μ$	Natural mortality rate
$β$	Disease transmission rate
$θ$	Rate at which exposed individuals become infectious
$δ$	Ebola-induced mortality rate
$c$	Proportion of cases that recover in isolation
$τ$	Total recovery rate (sum of isolated and non-isolated recovery)
$ρ$	Rate at which non-isolated recovered individuals transition back to infectious status due to viral persistence, enabling potential transmission to susceptible individuals

Figure 2. Time series (with and without delay) for tau2 = 1; lambda = 1; mu = 0.1; beta = 0.5; mu2 = 0.1; m = 0.5; p = 0.05; delta = 0.2; c = 0.7; rho = 0.1; tau = 0.6; theta = 0.3; S0 = 10; I0 = 0.1; R10 = 0; R20 = 0; tspan = [0, 200].

Figure 3. Time series (with and without delay) for tau2 = 2; lambda = 1; mu = 0.1; beta = 0.5; mu2 = 0.1; m = 0.5; p = 0.05; delta = 0.2; c = 0.7; rho = 0.1; tau = 0.6; theta = 0.3; S0 = 10; I0 = 0.1; R10 = 0; R20 = 0; tspan = [0, 200].

Figure 4. Time series (with and without delay) for tau2 = 3; lambda = 1; mu = 0.1; beta = 0.5; mu2 = 0.1; m = 0.5; p = 0.05; delta = 0.2; c = 0.7; rho = 0.1; tau = 0.6; theta = 0.3; S0 = 10; I0 = 0.1; R10 = 0; R20 = 0; tspan = [0, 200].

Figure 5. Time series (with and without delay) for tau2 = 4; lambda = 1; mu = 0.1; beta = 0.5; mu2 = 0.1; m = 0.5; p = 0.05; delta = 0.2; c = 0.7; rho = 0.1; tau = 0.6; theta = 0.3; S0 = 10; I0 = 0.1; R10 = 0; R20 = 0; tspan = [0, 200].

Figure 6. Time series (with and without delay) for tau2 = 5; lambda = 1; mu = 0.1; beta = 0.5; mu2 = 0.1; m = 0.5; p = 0.05; delta = 0.2; c = 0.7; rho = 0.1; tau = 0.6; theta = 0.3; S0 = 10; I0 = 0.1; R10 = 0; R20 = 0; tspan = [0, 200].

Figure 7. Time series (with and without delay) for tau2 = 10; lambda = 1; mu = 0.1; beta = 0.5; mu2 = 0.1; m = 0.5; p = 0.05; delta = 0.2; c = 0.7; rho = 0.1; tau = 0.6; theta = 0.3; S0 = 10; I0 = 0.1; R10 = 0; R20 = 0; tspan = [0, 200].

Figure 8. Time series (with and without delay) for both R_1 and R_2 where tau2 = 3; lambda = 1; mu = 0.1; beta = 0.5; mu2 = 0.1; m = 0.5; p = 0.05; delta = 0.2; c = 0.7; rho = 0.1; tau = 0.6; theta = 0.3; S0 = 10; I0 = 0.1; R10 = 0; R20 = 0; tspan = [0, 200].

Figure 9. Relationship between basic reproduction number R_0 and the isolation parameter c for R_0 > 1.

Figure 10. Relationship between basic reproduction number R_0 and the isolation parameter c for R_0 < 1.

Figure 11. How each parameter contribute to the basic reproduction number.

Figure 12. Phase space of equilibrium E* (S, , R*) when the delay is 1.5 i.e tau2 = 1.5. for lambda = 1; mu = 0.1; beta = 0.5; mu2 = 0.1; m = 0.05; delta = 0.2; c = 0.7; rho = 0.1; p = 0.05; S0 = 10; I0 = 0.1; R10 = 0; R20 = 0; tspan = [0, 200].

Figure 13. EE = [2.2313 1.3254 5.5867], DFE = [13 0 0]. Lambda = 1.30; mu = 0.1; beta = 0.5; delta = 0.2; rho = 0.1; c = 0.7; tau = 0.6; tau2 = 2.5; m = 0.05; mu2 = 0.1; S0 = 5; I0 = 1; R10 = 0; R20 = 0; tspan = [0, 300].

Figure 14. EE = [2.2313 1.3254 1.1963], DFE = [13 0 0]. Lambda = 1.30; mu = 0.1; beta = 0.5; delta = 0.2; rho = 0.1; c = 0.7; tau = 0.6; tau2 = 2.5; m = 0.05; mu2 = 0.1; S0 = 5; I0 = 1; R10 = 0; R20 = 0; tspan = [0, 300].

Figure 15. Lambda = 1.30; mu = 0.1; beta = 0.5; delta = 0.2; rho = 0.1; c = 0.7; tau = 0.6; tau2 = 2.5; m = 0.05; mu2 = 0.1; S0 = 5; I0 = 1; R10 = 0; R20 = 0; tspan = [0, 300].

Figure 16. Lambda = 1.30; mu = 0.1; beta = 0.5; delta = 0.2; rho = 0.1; c = 0.7; tau = 0.6; tau2 = 2.5; m = 0.05; mu2 = 0.1; S0 = 5; I0 = 1; R10 = 0; R20 = 0; tspan = [0, 300].

Figure 17. Lambda = 1.30; mu = 0.1; beta = 0.5; delta = 0.2; rho = 0.1; c = 0.7; tau = 0.6; tau2 = 2.5; m = 0.05; mu2 = 0.1; S0 = 5; I0 = 1; R10 = 0; R20 = 0; tspan = [0, 300].

8. Conclusion

This study provides a thorough analysis of the Ebola transmission dynamics using a deterministic SIRR model that incorporates key features such as nonlinear incidence rates and delayed recovery. By examining the stability of disease equilibria and the role of the basic reproduction number ( $R_{0}$ ), we have highlighted critical factors that influence the onset and persistence of outbreaks. The sensitivity analysis and bifurcation analysis, along with numerical simulations, underscore the importance of timely interventions and targeted strategies to reduce transmission, limit susceptible populations, and improve recovery rates. These insights offer valuable guidance for public health policy, particularly in managing outbreaks by identifying threshold conditions and critical parameters that determine the trajectory of disease spread. Overall, our findings contribute to a deeper understanding of Ebola dynamics and demonstrate the importance of considering delays in disease models to capture the complexities of real-world transmission.

Ethics and consent

This study did not involve human participants or personal data; therefore, ethical approval and consent were not required.

Data availability statement

As outlined in our manuscript, the present study is purely theoretical in scope, focusing on conceptual development and does not involve the generation or analysis of empirical data.

References

1. Kermack WO, McKendrick AG: A contribution to the mathematical theory of epidemics. Proc. R. Soc. Lond. A. 1927; 115(772): 700–721.
2. Brauer F: Mathematical Models in Epidemiology. Springer; 2017.
3. Anderson RM, May RM: The Population Dynamics of Infectious Diseases: Theory and Applications. Springer; 1981.
4. Erinle-Ibrahim LM, Idowu KO: Mathematical modelling of pneumonia dynamics of children under the age of five. Mathematical Association of Nigeria (MAN). 2021; 48(2): 119.
5. Hethcote HW: The mathematics of infectious diseases. SIAM Rev. 2000; 42(4): 599–653. Publisher Full Text
6. Heesterbeek H, Dietz K: The concept of RO in epidemic theory. Stat. Methods Med. Res. 1996; 5(4): 361–387.
7. Lekone P, Finkenstädt B: Statistical inference in a stochastic epidemic SEIR model with control intervention: Ebola as a case study. Biometrics. 2006; 62: 1170–1177. PubMed Abstract | Publisher Full Text
8. Ndanguza D, Tchuenche J, Haario H: Statistical data analysis of the 1995 Ebola outbreak in the Democratic Republic of Congo. Afr. Mat. 2013; 24: 55–68. Publisher Full Text
9. Fasina FO, et al.: Transmission dynamics and control of Ebola virus disease outbreak in Nigeria, July to September 2014. Euro Surveill. 2014; 19(40): 20920. PubMed Abstract | Publisher Full Text
10. Agusto FB, Teboh-Ewungkem MI, Gumel AB: Mathematical assessment of the effect of traditional beliefs and customs on the transmission dynamics of the 2014 Ebola outbreaks. BMC Med. 2015; 13: 1–17. Publisher Full Text
11. Idowu KO, Loyinmi AC: Impact of contaminated surfaces on the transmission dynamics of corona virus disease (Covid-19).2023.
12. Weitz JS, Dushoff J: Modeling post-death transmission of Ebola: challenges for inference and opportunities for control. Sci. Rep. 2015; 5: 1–7. Publisher Full Text
13. Ivorra B, Ngom D, Ramos Á: Be-CoDiS: A mathematical model to predict the risk of human diseases spread between countries—validation and application to the 2014–2015 Ebola virus disease epidemic. Bull. Math. Biol. 2015; 77: 1668–1704. PubMed Abstract | Publisher Full Text
14. Berge T, Lubuma JMS, Moremedi GM, et al.: A simple mathematical model for Ebola in Africa. J. Biol. Dyn. 2017; 11: 42–74. PubMed Abstract | Publisher Full Text
15. Rachah A, Torres DFM: Dynamics and optimal control of Ebola transmission. Math. Comput. Sci. 2016; 10: 331–342. Publisher Full Text
16. Abo SMC, Smith R: Modelling the daily risk of Ebola in the presence and absence of a potential vaccine. Infect. Dis. Model. 2020; 5: 905–917. PubMed Abstract | Publisher Full Text | Free Full Text
17. Marais F, et al.: A community-engaged infection prevention and control approach to Ebola. Health Promot. Int. 2016; 31: 440–449. PubMed Abstract | Publisher Full Text
18. Leroy EL, et al.: Fruit bats as reservoirs of Ebola virus. Nature. 2005; 438: 575–576. PubMed Abstract | Publisher Full Text
19. Burkett-Cadena ND, Vittor AY: Deforestation and vector-borne disease: Forest conversion favors important mosquito vectors of human pathogens. Sci. Rep. 2021; 11(1): 1–10.
20. Kucharski AJ, et al.: Measuring the impact of Ebola control measures in Sierra Leone. Proc. Natl. Acad. Sci. 2015; 112(46): 14366–14371. PubMed Abstract | Publisher Full Text | Free Full Text
21. Nishiura H, Chowell G: Transmission dynamics and control of Ebola virus disease: a review. BMC Med. 2014; 12(1): 196. Publisher Full Text
22. Agbomola JO, Lasekan AEEpidemiological model of COVID-19 with quarantine using spatially extended epidemiological model. Franklin Open. 2025; 12: 100327.PubMed Abstract | Publisher Full Text | Free Full Text
23. Oluwatobi IK, Chris LA: Qualitative analysis of the transmission dynamics and optimal control of covid-19.2023.
24. Legrand J, et al.: Understanding the dynamics of Ebola epidemics. Epidemiol. Infect. 2007; 135(4): 610–621. PubMed Abstract | Publisher Full Text | Free Full Text
25. Pandey A, et al.: Strategies for containing Ebola in West Africa. Science. 2014; 346(6212): 991–995. PubMed Abstract | Publisher Full Text | Free Full Text
26. Rivers CM, et al.: Modeling the impact of interventions on an epidemic of Ebola in Sierra Leone and Liberia. PLoS Curr. 2014; 6. Publisher Full Text
27. Althaus CL: Estimating the reproduction number of Ebola virus (EBOV) during the 2014 outbreak in West Africa. PLoS Curr. 2014; 6. PubMed Abstract | Publisher Full Text | Free Full Text
28. Loyinmi AC, et al.: On the interaction of the human immune system with foreign body: mathematical modeling approach. Kathmandu University Journal of Science, Engineering and Technology. 2023; 17(2). Publisher Full Text
29. Chretien J-P, et al.: Drought-associated chikungunya emergence along coastal East Africa. Am. J. Trop. Med. Hyg. 2007; 76(3): 405–407. PubMed Abstract | Publisher Full Text
30. Cooke KL, van den Driessche P : Analysis of an SEIRS epidemic model with two delays. J. Math. Biol. 1996; 35(2): 240–260. PubMed Abstract | Publisher Full Text
31. Kuang Y: Delay Differential Equations with Applications in Population Dynamics. Academic Press; 1993.
32. Smith HL: An Introduction to Delay Differential Equations with Applications to the Life Sciences. Springer; 2011.
33. Hale JK, Lunel SMV: Introduction to Functional Differential Equations. Springer-Verlag; 1993.
34. Beretta E, Kuang Y: Geometric stability switch criteria in delay differential systems with delay-dependent parameters. SIAM J. Math. Anal. 2002; 33(5): 1144–1165. Publisher Full Text
35. Shampine LF, Thompson S: Solving DDEs in MATLAB. Appl. Numer. Math. 2000; 37(4): 441–458. Publisher Full Text
36. Olawale-Shosanya SO, et al.: A Meta-Ensemble Predictive Model For The Risk Of Lung Cancer. Al-Bahir J. Eng. Pure Sci. 2024; 5(1): 4. Publisher Full Text
37. Idowu KO, et al.: A mathematical model investigating the impact of campaign program on the pattern of transmission of anemia in pregnant women. Kathmandu University Journal of Science, Engineering and Technology. 2024; 18(1): 1–12. Publisher Full Text
38. Erinle-Ibrahim LM, et al.: Application of homotopy perturbation method to the mathematical modelling of temperature rise during microwave hyperthermia. Fudma J. Sci. 2021; 5(2): 273–282. Publisher Full Text
39. Agbomola JO, Loyinmi AC: Modelling the impact of some control strategies on the transmission dynamics of Ebola virus in human-bat population: An optimal control analysis. Heliyon. 2022; 8(12): e12121. PubMed Abstract | Publisher Full Text | Free Full Text
40. Castillo-Chavez C, Song B: Dynamical models of tuberculosis and their applications. Math. Biosci. Eng. 2004; 1(2): 361–404. Publisher Full Text
41. Huang G, Takeuchi Y, Ma W: Lyapunov functionals for delay differential equations model of viral infections. SIAM J. Appl. Math. 2010; 70(7): 2693–2708. Publisher Full Text
42. Sastry S: Nonlinear Systems: Analysis, Stability, and Control. Springer; 1999.

Comments on this article Comments (0)

Version 1

VERSION 1 PUBLISHED 02 Sep 2025