Kelachi Enwere^*, Uchenna Ogoke, Ethelbert Nduka
Department of Mathematics and Statistics, University of Port Harcourt, Port Harcourt, Nigeria.
DOI: 10.4236/oalib.1114397   PDF    HTML   XML   46 Downloads   458 Views   Citations

Abstract

Infectious diseases such as bacterial, viral, fungal and parasitic diseases continue to have a significant impact on global health. Hence, knowledge on how infection spread, population vulnerable to it and factors that determine where and when these infectious diseases occur is crucial for improving clinical trials and proffer effective treatment. This research work seeks to model disease incidence of malaria using the SEIR-SEI compartmental model. Secondary data on the incidence of Malaria disease from the World Malaria Report for Nigeria for 2023, Coefficients from Macro trends, were used to study the behavior of the Compartmental SEIR-SEI model. The results from the study show that since R < 1, then disease-free equilibrium for malaria transmission can be achieved as the number of infected individuals eventually declines in the population. Our herd immunity threshold decreases to about 0.3582 at t = 4 which means that the number of infected individuals decreases as a result of increase in vaccination of the population from the onset of disease incidence. Nevertheless, there will always be new infectious viruses, bacteria, fungal and parasites with their corresponding outbreaks. However, Compartmental modelling on how infectious diseases progress in an epidemic can inform how effective intervention will be.

Keywords

Susceptible, Infectious, Malaria, Exposed, Compartment, Immunity

Share and Cite:

1. Introduction

Throughout recorded history, human society has lived with periodic epidemics and pandemics. Numerous disease outbreaks have led to death, societal upheaval, and economic disruption. Predicting how an outbreak may progress is therefore essential to mitigate its effects, and the field of epidemiologic modeling is central to this.

Infectious diseases progress within populations both due to the behavior of the infectious agent and the population itself. Models of how they progress in an epidemic are based on a set of assumptions and statistics, which are used to establish a set of parameters that inform how effective intervention will be (for example, social distancing or mass vaccination). This can be used to predict which interventions to implement or avoid as well as future growth and spread patterns and many other variables.

Infectious diseases such as bacterial, viral, fungal and parasitic diseases continue to have a significant impact on global health. Many new infectious diseases such as HIV, Ebola virus, and the coronavirus disease of 2019 (COVID-19) have also emerged to pose additional health threats at global and national levels. Hence, a proper understanding of disease dynamics is crucial for monitoring a vaccination program and assessing the risk of future outbreaks. This might be particularly important in resource-constrained settings and in countries with a high infectious diseases burden.

Consequently, information learning techniques and models such as compartmental modelling have been developed to aid those in the field of medicine, epidemiology and biostatistics to investigate the effects infectious diseases have on the population of interest. It has also been used to identify potential hazards, exposure and vulnerabilities in order to improve clinical trials and proffer effective treatment. In addition, such models can help interpret data analysis on infectious disease and improve control measures.

Nevertheless, knowledge on how infection spread, population vulnerable to it and factors that determine where and when these infectious diseases occur is crucial for policy makers to implement targeted interventions in order to mitigate its effects. See Figure 1 below on the timeline for an infection and its adverse effects on individuals of the population.

Figure 1. Natural history timeline for an infection.

Figure 1 shows the dynamics of infectiousness for a contagious disease and its effect of spread at various stages. For at an initially detection of an infection, a proportion of the population becomes susceptible to the disease due to random mixing, which could spread over time. Individuals who may have been exposed to the disease may show signs of disease symptoms and are eventually infected. This proportion of infectious individuals, sick of the disease may then recover through treatments or vaccination if they escape mortality. Moreover, the results could be death from the infection or recovery from treatments or vaccination and/or are susceptible again to the infection.

1.1. Malaria

Malaria is a life-threatening disease spread to humans through the bites of an infectious female Anopheles mosquito. It is mostly found in tropical countries with prevailing environmental condition of warm humid climate, swampy environments, poor sanitary habits, poverty and ignorance exist.

Globally in 2023, there were an estimated 263 million malaria cases and 597,000 malaria deaths in 83 countries. In 2023, the World Health Organization (WHO) stated that “African Region was home to 94% of malaria cases (246 million) and 95% (569,000) of malaria deaths. Children under 5 accounted for about 76% of all malaria deaths in the Region. Over half of these deaths occurred in four countries: Nigeria (30.9%), the Democratic Republic of the Congo (11.3%), Niger (5.9%) and United Republic of Tanzania (4.3%)” [1].

1.2. Compartmental Model

Compartmental models are often applied to the mathematical modelling of infectious diseases. These models are used to analyze the disease dynamics and to estimate the total number of infected people, the total number of recovered people, and to estimate epidemiological parameters such as the basic reproduction number or effective reproduction number. Such models can show how different public health interventions may affect the outcome of an epidemic. The models are conceptually simple and thus very useful in the field of epidemiology, medicine, economics, genetics, engineering, mathematics etc.

A compartment model partition the population into different compartments (or classes) such as Susceptible (S), Exposed (E), Infected (I), Hospitalized (H), Quarantined (Q), Recovered (R), and Death (D) compartments. The population is assigned to the compartments and people may progress between compartments with change over time t. The acronym for any compartmental model usually show how individuals moves within the compartments; for example, SEIS means susceptible, exposed, infectious, then susceptible again. Compartmental models are most often run with ordinary differential equations (which are deterministic), but can also be used with a stochastic (random) framework, which is more realistic but much more complicated to analyze.

One of the most common compartmental models is the Susceptibility, Infection, and Recovery Models (SIR). The SIR model accurately represents how an infection would spread through a population because it takes into consideration that some people will recover from the disease and no longer be susceptible. Hence, this model assumes that people who recover from the infection become immune and cannot become infected a second time.

Compartmental models may be used to predict properties of how a disease spreads, for example the prevalence (total number of infected) or the duration of an epidemic. In addition, the model allows for understanding how different situations may affect the outcome of the epidemic, e.g., what the most efficient technique is for issuing a limited number of vaccines in a given population.

1.3. Statement of the Problem

Throughout recorded history, human society has lived with periodic epidemics and pandemics. Infectious diseases such as bacterial, viral, and parasitic diseases continue to have a significant impact on global health. Numerous disease outbreaks have led to death, societal upheaval, and economic disruption. Hence, a proper understanding of disease dynamics is crucial for monitoring a vaccination program and assessing the risk of future outbreaks.

Nevertheless, in recent years, there has been an increase in research activity regarding compartmental modelling for a complete understanding of disease transmission and persistence. This study seeks to determine the dynamics of spread on the incidence of malaria disease, its infectiousness and control measures using the SEIR-SEI compartmental model.

1.4. Overview of Literature

Mathematical models have a long history of being used to describe the spread of infectious diseases from plague outbreaks more than a century ago to the more recent SARS and Ebola epidemics [2]. From making decisions around different vaccination strategies for influenza to modelling HIV, and from modelling pandemic influenza to currently facilitating real-time policy decision making around the COVID-19 epidemic. However, there is increasing research on mathematical models to fully understand the transmission and persistence of infectious diseases. With ongoing advances in computational tools as well as access to disease incidence data, the use of such models continues to increase.

[3] as well as [4] studied the dynamical behaviour of epidemiological models with nonlinear incidence rates. Their studies show that models with nonlinear incidence have a much wider range of dynamical behaviours than do those with bilinear incidence rates. Again, [5] studied the global dynamics of a SEIR model with varying total population size.

[6] investigated the impact of immigration on the transmission dynamics of tuberculosis. Their theoretical analysis indicated that the disease would persist in the population if there were a prevalence of TB in immigrants. Their study suggests that immigrants have a considerable influence on the overall transmission dynamics behaviour of tuberculosis. Their study show that when there is an immigration of infected people into a population, then the disease-free equilibrium cannot be achieve. [7] incorporated treatment of individuals infectious to TB and chemoprophylaxis (treatment for the latently infected) on the SEIR model. The model assumed that the latently infected individuals develop active disease because of endogenous re-activation, exogenous re-infection and disease relapse. The study shows that treatment of infectious individuals is more effective in the first years of implementation as it cleared active TB immediately. As a result, chemoprophylaxis will do better in controlling the number of infectious due to reduced progression to active TB.

Nevertheless, [8] used the SIR model to describe heterogeneity in vaccination and the cost effectiveness of supplemental immunization activities for measles in Uganda. [9] also developed the SIR model to simulate and better understand the multi-periodic patterns in outbreaks of avian flu in North America. In addition, [10] used a compartmental SEIR model to simulate the dynamics of Ebola outbreak in the Democratic Republic of Congo in 1995. [11] proposed a stochastic discrete time variable of the SIR model for infectious diseases transmission within and between districts, and susceptible individuals interact.

However, Public health organizations throughout the world use such models to evaluate and develop intervention disease outbreak policies for ever-emerging epidemics. Simulation allows for rapid assessment and decision-making, providing quantification and insight into the dynamics of a spread. An intensive inter- and multidisciplinary research effort is speeding up the developments in the field integrating advances from epidemiology, molecular biology, computational engineering and science, and applied mathematics as well as sociology.

To this end, this study seeks to model the dynamics of an infectious disease, using the SEIR-SEI compartmental Model. In addition, to illustrate disease infectiousness to control measures using model on Malaria as case study.

2. Methodology

For the purpose of this research, secondary data on the incidence of Malaria disease from the “World Malaria Report for 2023, by the World Health Organization for Nigeria” [1], Coefficients from Macro trends [12], were used to study the behavior of the Compartmental SEIR-SEI model.

2.1. Susceptibility, Exposed, Infection, Recovery and Susceptible, Exposed, Infected (SEIR-SEI) Compartmental Model

SEIR-SEI compartmental model describes interaction between susceptible, exposed, infected and recovered human population and the susceptible, exposed and infected mosquito population (Figure 2). The SEIR-SEI model is evaluated using partial differential equations (which are deterministic), but can also be used with a stochastic (random) framework, which is more realistic but much more complicated to analyze. Individuals in the human population as well as mosquitoes, may move between compartments with change over time t. The SEIR human model part can be used when vital dynamics (births and deaths) are enabled which sustain an epidemic or allow new introductions to spread. Transmission rate depends on the prevalence of infective, the contact rate (biting rate) and the probability of transmission given contact.

Figure 2. Diagram of a SEIR-SEI compartmental model with birth and death (for human and mosquito population).

Table 1 below defines the above model parameters as well as their numerical value and source.

Table 1. Parameter description for SEIR-SEI compartmental model with birth and death (for human and mosquito population).

PARAMETER	DESCRIPTION	NUMERICAL VALUE	SOURCE
${\lambda }_h$	Natural birth rate in host population	0.035341	Micro trends
$S_h$	Susceptible individuals amongst the population who may be randomly beaten by an infected mosquito.	2185	WHO
${\mu }_h$	Natural death rate for host population.	0.010672	Micro trends
${\beta }_h$	Contact rate of infection from vector to host.	0.00035	Assumed
$E_h$	Exposed individuals to mosquito bites but not yet infectious.	1741	WHO
$E_R$	Exposed individuals safe from being infected due to natural immunity or through vaccination.
${\sigma }_h$	Progression rate from being exposed to being infected.	0.00031	Assumed
$I_h$	Infected individuals from malaria.	681	WHO
${\mu }_I$	Mortality from death of infection.
${\gamma }_h$	Recovery rate in host population.	0.071429	Calculated
$R_h$	Individuals who recover from the infection in the host population.	1057	WHO
$R_S$	Individuals in the host population, after recovery are susceptible again to the disease (via mosquito bites).
${\lambda }_m$	Birth rate in the vector population.		Not Considered
$S_m$	Number of susceptible in the vector population that may be infected with malaria disease.		Not Considered
${\mu }_m$	Death rate in the vector population.		Not Considered
${\beta }_m$	Contact rate of infection from host to vector.		Not Considered
$E_m$	Exposed mosquito to an infected human.		Not Considered
${\sigma }_m$	Progression rate of vector from being exposed to being infected.		Not Considered
$I_m$	Number of infected mosquitoes from malaria.		Not Considered

Hence, the SEIR-SEI model for malaria disease for the population (host and vector) of $N_h$ individuals and $N_m$ mosquitoes consist of a set of seven differential equations where the solution functions $S_h$ , $E_h$ , $I_h$ and $R_h$ for host and $S_m$ , $E_m$ and $I_m$ for the vector are single variable real-valued functions that $S_h$ , $E_h$ , $I_h$ , $R_h$ , $S_m$ , $E_m$ , $I_m$ : $R_0^{+}\to \left[0,1\right]$ and ${\beta }_h$ , ${\lambda }_h$ , ${\sigma }_h$ , ${\gamma }_h$ , ${\mu }_h$ , ${\beta }_m$ , ${\lambda }_m$ , ${\sigma }_m$ , ${\gamma }_m$ , ${\mu }_m\in R_0^{+}$ .

In addition, the total population proportion at any time t, is assumed to satisfy $S_h\left(t\right)+E_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)=1$ ; and $S_m\left(t\right)+E_m\left(t\right)+I_m\left(t\right)=1$

Thus, the SEIR-SEI model for the host-vector population are given as,

Susceptible, ($S_h$ ) ⇒ $\frac{\alpha S_h}{\alpha t}={\lambda }_h-{\beta }_h{S}_h{I}_h-{\mu }_h{S}_h$ (1)

Exposed, ($E_h$ ) ⇒ $\frac{\alpha E_h}{\alpha t}={\beta }_h{S}_h{I}_h-{\sigma }_h{E}_h-{\mu }_h{E}_h$ (2)

Infectious, ($I_h$ ) ⇒ $\frac{\alpha I_h}{\alpha t}={\sigma }_h{E}_h-{\gamma }_h{I}_h-{\mu }_h{I}_h$ (3)

Recovered, ($R_h$ ) ⇒ $\frac{\alpha R_h}{\alpha t}={\gamma }_h{I}_h-{\mu }_h{R}_h$ (4)

Susceptible, ($S_m$ ) ⇒ $\frac{\alpha S_m}{\alpha t}={\lambda }_m-{\beta }_m{S}_m{I}_m-{\mu }_m{S}_m$ (5)

Exposed, ($E_m$ ) ⇒ $\frac{\alpha E_m}{\alpha t}={\beta }_m{S}_m{I}_m-{\sigma }_m{E}_m-{\mu }_m{E}_m$ (6)

Infectious, ($I_m$ ) ⇒ $\frac{\alpha I_m}{\alpha t}={\sigma }_m{E}_m-{\gamma }_m{I}_m-{\mu }_m{I}_m$ (7)

where $\beta =\tau \overline{c}$ is the effective contact rate, and $\gamma$ is the removal rate. This means that the expected duration of infection is simply the inverse of the removal rate, i.e. $d={\gamma }^{-1}$

In this study, we are more interested in malaria transmission from vector to human. Therefore, from the above differential equations, we would be concerned on the susceptible, exposed infectious and recovered equations of the host population of human. Thus, we can solve the SEIR model part of the human population for a small change in time t, (or short step in time) using the Taylor method, Euler method and the Runge-Kutta method [13]-[15].

Consequently, the Taylor’s method becomes exceedingly complex as its order increases because each increase in order requires a further derivative to be taken for each of our differential equations in the system.

Recall that the Taylor’s expansion of order “n” is given as,

$y_{i+1}=y_i+hf\left(t_i,y_i\right)+\frac{h^2}{2!}{f}^{\prime }\left(t_i,y_i\right)+\frac{h^3}{3!}{f}^{″}\left(t_i,y_i\right)+\cdots +\frac{h^n}{n!}{f}^{\left(n-1\right)}\left(t_i,y_i\right)$ (8)

$i=0,1,2,\cdots ,n$ $h=\frac{t_{\max }-t_{\min }}{n}$ $t_{i+1}=t_{\min }+ih$ $t_0=t_{\min }$ $f=\frac{\text{d}y}{\text{d}t}$

where the dependent variable $y_i$ is calculated by short steps of equal time interval $h$ of the independent variable $x$ .

Hence, the Euler method is the Taylor method for order one, i.e.

$y_{i+1}=y_i+hf\left(t_i,y_i\right)$ (9)

Thus, the SEIR model with birth and death for the host population with dependent variables $S_h\left(t\right)$ , $E_h\left(t\right)$ , $I_h\left(t\right)$ and $R_h\left(t\right)$ are modeled using the Euler’s iterative method represented by the equations

$S_{hi+1}=S_{hi}+h\frac{\text{d}{S}_h}{\text{d}t}\left({\lambda }_h-{\beta }_h{S}_h{I}_h-{\mu }_h{S}_h\right)$ (10)

$E_{hi+1}=E_{hi}+h\frac{\text{d}{E}_h}{\text{d}t}\left({\beta }_h{S}_h{I}_h-{\sigma }_h{E}_h-{\mu }_h{E}_h\right)$ (11)

$I_{hi+1}=I_{hi}+h\frac{\text{d}{I}_h}{\text{d}t}\left({\sigma }_h{E}_h-{\gamma }_h{I}_h-{\mu }_h{I}_h\right)$ (12)

$R_{hi+1}=R_{hi}+h\frac{\text{d}{R}_h}{\text{d}t}\left({\gamma }_h{I}_h-{\mu }_h{R}_h\right)$ (13)

On the other hand, the Runge-Kutta method is the Taylor method of order four having four sections ($K_1$ , $K_2$ , $K_3$ , $K_4$ ) of equations. Thus, for an order four Taylor method, i.e.

$y_{i+1}=y_i+hf\left(t_i,y_i\right)+\frac{h^2}{2!}{f}^{\prime }\left(t_i,y_i\right)+\frac{h^3}{3!}{f}^{″}\left(t_i,y_i\right)$

the SEIR model with birth and death for the host population with dependent variables $S_h\left(t\right)$ , $E_h\left(t\right)$ , $I_h\left(t\right)$ and $R_h\left(t\right)$ are modeled using the Runge-Kutta’s iterative method represented by the equation

$y_{i+1}=y_i+\frac{1}{6}\left(k_1+2k_2+2k_3+k_4\right)$ (14)

where

$K_1=hf\left(t_i,y_i\right)$ $K_2=hf\left(t_i+\frac{h}{2},y_i+\frac{K_1}{2}\right)$ $K_3=hf\left(t_i+\frac{h}{2},y_i+\frac{K_2}{2}\right)$

$K_4=hf\left(t_{i+1},y_i+k_3\right)$ $t_{i+1}=t_{\min }+ih$

Since the SEIR model for host population is a system of four equations, the model will include sixteen separate equations given by,

$\begin{array}{l}{K}_{1S_h}=h\frac{\text{d}{S}_h}{\text{d}t}\left(t_i,S_i,I_i\right)\\ K_{1S_h}={\lambda }_h-{\beta }_h{S}_h{I}_h-{\mu }_h{S}_h\end{array}$ (15)

$\begin{array}{l}{K}_{1E_h}=h\frac{\text{d}{E}_h}{\text{d}t}\left(t_i,S_i,E_i,I_i\right)\\ K_{1E_h}={\beta }_h{S}_h{I}_h-{\sigma }_h{E}_h-{\mu }_h{E}_h\end{array}$ (16)

$\begin{array}{l}{K}_{1I_h}=h\frac{\text{d}{I}_h}{\text{d}t}\left(t_i,E_i,I_i\right)\\ K_{1I_h}={\sigma }_h{E}_h-{\gamma }_h{I}_h-{\mu }_h{I}_h\end{array}$ (17)

$\begin{array}{l}{K}_{1R_h}=h\frac{\text{d}{R}_h}{\text{d}t}\left(t_i,I_i\right)\\ K_{1R_h}={\gamma }_h{I}_h-{\mu }_h{R}_h\end{array}$ (18)

$\begin{array}{l}{K}_{2S_h}=h\frac{\text{d}{S}_h}{\text{d}t}\left(t_i+\frac{h}{2},S_{hi}+\frac{1}{2}{K}_{1S_{hi}},I_{hi}+\frac{1}{2}{K}_{1I_{hi}}\right)\\ K_{2S_h}={\lambda }_h-{\beta }_h\left(S_{hi}+\frac{1}{2}{K}_{1S_{hi}}\right)\left(I_{hi}+\frac{1}{2}{K}_{1I_{hi}}\right)-{\mu }_h\left(S_{hi}+\frac{1}{2}{K}_{1S_{hi}}\right)\end{array}$ (19)

$\begin{array}{l}{K}_{2E_h}=h\frac{\text{d}{E}_h}{\text{d}t}\left(t_i+\frac{h}{2},S_{hi}+\frac{1}{2}{K}_{1S_{hi}},E_{hi}+\frac{1}{2}{K}_{1E_{hi}},I_{hi}+\frac{1}{2}{K}_{1I_{hi}}\right)\\ K_{2E_h}={\beta }_h\left(S_{hi}+\frac{1}{2}{K}_{1S_{hi}}\right)\left(I_{hi}+\frac{1}{2}{K}_{1I_{hi}}\right)-{\sigma }_h\left(E_{hi}+\frac{1}{2}{K}_{1E_{hi}}\right)\\ -{\mu }_h\left(E_{hi}+\frac{1}{2}{K}_{1E_{hi}}\right)\end{array}$ (20)

$\begin{array}{l}{K}_{2I_h}=h\frac{\text{d}{I}_h}{\text{d}t}\left(t_i+\frac{h}{2},E_{hi}+\frac{1}{2}{K}_{1E_{hi}},I_{hi}+\frac{1}{2}{K}_{1I_{hi}}\right)\\ K_{2I_h}={\sigma }_h\left(E_{hi}+\frac{1}{2}{K}_{1E_{hi}}\right)-{\gamma }_h\left(I_{hi}+\frac{1}{2}{K}_{1I_{hi}}\right)-{\mu }_h\left(I_{hi}+\frac{1}{2}{K}_{1I_{hi}}\right)\end{array}$ (21)

$\begin{array}{l}{K}_{2R_h}=h\frac{\text{d}{I}_h}{\text{d}t}\left(t_i+\frac{h}{2},I_i+\frac{1}{2}{K}_{1I_{hi}}\right)\\ K_{2R_h}={\gamma }_h\left(I_i+\frac{1}{2}{K}_{1I_{hi}}\right)-{\mu }_h\left(R_i+\frac{1}{2}{K}_{1R_{hi}}\right)\end{array}$ (22)

$\begin{array}{l}{K}_{3S_h}=h\frac{\text{d}{S}_h}{\text{d}t}\left(t_i+\frac{h}{2},S_{hi}+\frac{1}{2}{K}_{2S_{hi}},I_{hi}+\frac{1}{2}{K}_{2I_{hi}}\right)\\ K_{3S_h}={\lambda }_h-{\beta }_h\left(S_{hi}+\frac{1}{2}{K}_{2S_{hi}}\right)\left(I_{hi}+\frac{1}{2}{K}_{2I_{hi}}\right)-{\mu }_h\left(S_{hi}+\frac{1}{2}{K}_{2S_{hi}}\right)\end{array}$ (23)

$\begin{array}{l}{K}_{3E_h}=h\frac{\text{d}{E}_h}{\text{d}t}\left(t_i+\frac{h}{2},S_{hi}+\frac{1}{2}{K}_{2S_{hi}},E_{hi}+\frac{1}{2}{K}_{2E_{hi}},I_{hi}+\frac{1}{2}{K}_{2I_{hi}}\right)\\ K_{3E_h}={\beta }_h\left(S_{hi}+\frac{1}{2}{K}_{2S_{hi}}\right)\left(I_{hi}+\frac{1}{2}{K}_{2I_{hi}}\right)-{\sigma }_h\left(E_{hi}+\frac{1}{2}{K}_{2E_{hi}}\right)\\ -{\mu }_h\left(E_{hi}+\frac{1}{2}{K}_{2E_{hi}}\right)\end{array}$ (24)

$\begin{array}{l}{K}_{3I_h}=h\frac{\text{d}{I}_h}{\text{d}t}\left(t_i+\frac{h}{2},E_{hi}+\frac{1}{2}{K}_{2E_{hi}},I_{hi}+\frac{1}{2}{K}_{2I_{hi}}\right)\\ K_{3I_h}={\sigma }_h\left(E_{hi}+\frac{1}{2}{K}_{2E_{hi}}\right)-{\gamma }_h\left(I_{hi}+\frac{1}{2}{K}_{2I_{hi}}\right)-{\mu }_h\left(I_{hi}+\frac{1}{2}{K}_{2I_{hi}}\right)\end{array}$ (25)

$\begin{array}{l}{K}_{3R_h}=h\frac{\text{d}{I}_h}{\text{d}t}\left(t_i+\frac{h}{2},I_i+\frac{1}{2}{K}_{2I_{hi}}\right)\\ K_{3R_h}={\gamma }_h\left(I_{hi}+\frac{1}{2}{K}_{2I_{hi}}\right)-{\mu }_h\left(R_{hi}+\frac{1}{2}{K}_{2R_{hi}}\right)\end{array}$ (26)

$\begin{array}{l}{K}_{4S_h}=h\frac{\text{d}{S}_h}{\text{d}t}\left(t_{i+1},y_i+k_3\right)\\ K_{4S_h}={\lambda }_h-{\beta }_h\left(S_{hi}+K_{3S_{hi}}\right)\left(I_{hi}+K_{3I_{hi}}\right)-{\mu }_h\left(S_{hi}+K_{3S_{hi}}\right)\end{array}$ (27)

$\begin{array}{l}{K}_{4E_h}=h\frac{\text{d}{E}_h}{\text{d}t}\left(t_{i+1},y_i+k_3\right)\\ K_{4E_h}={\beta }_h\left(S_{hi}+K_{3S_{hi}}\right)\left(I_{hi}+K_{3I_{hi}}\right)-{\sigma }_h\left(E_{hi}+K_{3E_{hi}}\right)-{\mu }_h\left(E_{hi}+K_{3E_{hi}}\right)\end{array}$ (28)

$\begin{array}{l}{K}_{4I_h}=h\frac{\text{d}{I}_h}{\text{d}t}\left(t_{i+1},y_i+k_3\right)\\ K_{4I_h}={\sigma }_h\left(E_{hi}+K_{3E_{hi}}\right)-{\gamma }_h\left(I_{hi}+K_{3I_{hi}}\right)-{\mu }_h\left(I_{hi}+K_{3I_{hi}}\right)\end{array}$ (29)

$\begin{array}{l}{K}_{4R_h}=h\frac{\text{d}{R}_h}{\text{d}t}\left(t_{i+1},y_i+k_3\right)\\ K_{4R_h}={\gamma }_h\left(I_{hi}+K_{3I_{hi}}\right)-{\mu }_h\left(R_{hi}+K_{3R_{hi}}\right)\end{array}$ (30)

Thus, in general, Runge-Kutta method of order four SEIR model yields,

$S_{hi+1}=S_{hi}+\frac{1}{6}\left(K_{1S_{hi}}+2K_{2S_{hi}}+2K_{3S_{hi}}+K_{4S_{hi}}\right)$ (31)

$E_{hi+1}=E_{hi}+\frac{1}{6}\left(K_{1E_{hi}}+2K_{2E_{hi}}+2K_{3E_{hi}}+K_{4E_{hi}}\right)$ (32)

$I_{hi+1}=I_{hi}+\frac{1}{6}\left(K_{1I_{hi}}+2K_{2I_{hi}}+2K_{3I_{hi}}+K_{4I_{hi}}\right)$ (33)

$R_{hi+1}=R_{hi}+\frac{1}{6}\left(K_{1R_{hi}}+2K_{2R_{hi}}+2K_{3R_{hi}}+K_{4R_{hi}}\right)$ (34)

Furthermore, from both iterative methods, the results produced by the Runge-Kutta’s method is more accurate than that of the Euler’s method, since it has a forth order accuracy. [13]

2.2. Reproductive Number R

The Reproductive Number R is the average number of new infections that one infectious person produces. It measures how many people an infected individual will transmit the disease to before they recover. At the beginning of an epidemic when nearly all are susceptible, it is called R₀. Thus, the reproductive number R refers to the average number of people an infectious person will infect, assuming that the rest of the population is susceptible.

The reproduction number depends on biological, behavioural, social and environmental factors that could foster the spread of the infectious disease. When $R>1$ , an epidemic is most likely to occur. But if $R<1$ , there could still be (small) cluster of cases, but not enough to sustain an epidemic. In addition, $R=1$ is a threshold between epidemic and no epidemic. Note that for a given model, $R_0$ is fixed over all time.

Figure 3. Graphical depiction of Generations in an epidemic.

This R number is important in determining how quickly an infection will spread through a population. For example, if the value of R for measles is 14, then each case of measles would produce 14 new secondary cases. This would spread through the population much faster than where the value of R was only equal to 2.

Figure 3. describes generations in an epidemic, where the number of infections multiplies by a constant factor each day. Here, the reproductive number for the infectiousness above, $R=3$ , i.e. an infected person produces three new infections each day and this follows an exponential growth through each generations. This can be used as a simplistic model of how an infection could potentially spread through a population.

However, while the calculation of Reproductive Number R assumes that the entire population is susceptible to the disease, this may not always be the case. Since some people will be immune to the disease due to a prior infection creating life-long immunity, or because of vaccination. Hence, it is imperative to maintain knowledge of the value of Reproductive Number R, in order to control and potentially eradicate a disease from a given population. Thus, at each time period, the effective Reproduction Number $R_E$ is calculated as,

$R_E=Rs$ (35)

where $s$ is the fraction of the host population who are susceptible at time t.

Consequently, our R for SEIR model can be expressed as,

$R=\frac{\text{Infection}}{\text{Contact}}\times \frac{\text{Contact}}{\text{Time}}\times \frac{\text{Time}}{\text{Infection}}\times \frac{\text{Exposure}}{\text{Exposure}}$

i.e. $R=\tau \times \overline{c}\times d\times \epsilon$ (36)

where $\tau$ is the transmission probability of infection per contact between a susceptible and infected individual, $\overline{c}$ is the average rate of contact between susceptible and infected individuals, and $d$ is the duration of infectious period. $\epsilon$ is the probability of surviving exposure to infection and the inverse of being infected from the exposure, i.e. $\epsilon =\sigma {\sigma }^{-1}$

Thus, substituting for the effective contact rate $\beta =\tau \overline{c}$ and for the expected duration of infection $d={\gamma }^{-1}$ we have,

$R=\frac{\beta \epsilon }{\gamma \epsilon }$ (37)

Consequently, when $R<1$ then disease-free equilibrium can be achieved and when it is stable, the number of infected individuals eventually declines to zero, and the disease dies out from the population. On the other hand, the endemic equilibrium is asymptotically stable only when $R>1$ , meaning the infection persist in the population because each infected individual infects more than one other person on average.

2.3. Herd Immunity Threshold v

Herd immunity occurs when a critical proportion of a susceptible population is immunised against a contagious (infectious) disease, giving overall protection to the remainder of the “unprotected” herd (community). Herd immunity works because it is more difficult for diseases to spread between individuals if large numbers are already immune, as this breaks the chain of infectiousness.

Therefore, Herd Immunity Threshold $v$ , is the proportion of the population that needs to be vaccinated or immune in order for an infection to be contained. This is reached when each new case leads to just a single new case, leading to the infection becoming stable within the population. i.e. $R_E=1$ .

Consequently, at the onset of an infection, everyone in the population is either vaccinated or they are susceptible, so $s+v=1$ ; this implies $s=1-v$ and substituting this into (35), we have,

$R\left(1-v\right)=1$ (38)

Thus, Herd Immunity Threshold

$v=1-\frac{1}{R}$ (39)

where R is the reproduction number at each period.

2.4. Numerical Illustration on Malaria Using SEIR-SEI Model

The spread of Malaria disease is illustrated below from the given SEIR-SEI model part for the human population. The model employs data gotten from the “World Malaria Report for 2023, by the World Health Organization for Nigeria”. The incidence of Malaria is the estimated number of new and relapse Malaria cases arising in a given year, expressed as the rate per 100,000 population based on WHO World Malaria Report.

The SEIR model for human host, defined the Nigerian Population Size $N_h$ as 2279 the total population per 100,000 people. The susceptible group $S_h$ is the population of vulnerable people (2185 per 100,000 people) nonimmune to the disease. $E_h$ represents the exposed group of individuals, which means the individual has been beaten by mosquitoes but has not been infected malaria disease yet. Thus, the estimated number of individuals exposed to malaria in 2023 is (1741 per 100,000 people). The Infective $I_h$ is a group of individuals (681 per 100,000 people) that are sick of malaria disease which are to be treated, vaccinated (or perhaps quarantined or hospitalized so as not to spread the disease until they recover). $R_h$ is a group of recovered individuals that come from the infective group, having escaped death through treatment. WHO malaria report estimated 1057 per 100,000 people who recovered from malaria illness in 2023.

The population is open with Nigerian birth rate ${\lambda }_h$ of 0.035341 and mortality rate ${\mu }_h$ of 0.010672 respectively. Nevertheless, WHO duration of illness for an infected person diagnosed with malaria, to recovery may take up to about 2 weeks (14 days) into course of treatment. Thus, the removal rate ${\gamma }_h=\frac{1}{14}=0.071429$ . If the transmission rate of contact (biting rate) ${\beta }_h$ = 0.00035 and the progression rate ${\sigma }_h$ from exposures to mosquito bites is 0.00031, then, the SEIR model for the spread of malaria disease amongst human population can be modelled as follows:

As at the time of disease incidence (WHO 2023 world malaria report) $N_h=S_h+E_h+I_h+R_h$ (per 100,000 people). $S_h\left(0\right)=2185$ , $E_h\left(0\right)=1741$ , $I_h\left(0\right)=681$ and $R_h\left(0\right)=1057$ .

Therefore, we proceed to approximate the rates of change at the initial stage of identification at time t = 0 using the Runge-Kutta’s method.

3. Results and Discussion

For each time t, we proceed to approximate the rates of change at the initial stage of identification at time t = 0 using the fourth order Runge-Kutta iterative formula. Thus, modelling the SEIR human population compartments, we have

$K_{1S_h}=0.035341-\left[0.00035\times 2185\times 681\right]-\left[0.010672\times 2185\right]=-544\left(\text{persons}\right)$

$K_{1E_h}=\left[0.00035\times 2185\times 681\right]-\left[0.00031\times 1741\right]-\left[0.010672\times 1741\right]=502$

$K_{1I_h}=\left[0.00031\times 1741\right]-\left[0.071429\times 681\right]-\left[0.010672\times 681\right]=-55$

$K_{1R_h}=\left[0.071429\times 681\right]-\left[0.010672\times 1057\right]=37$

$\begin{array}{c}{K}_{2S_h}=0.035341-0.00035\left(2185+\frac{-544}{2}\right)\left(681+\frac{-55}{2}\right)-0.010672\left(2185+\frac{-544}{2}\right)\\ =0.035341-0.00035\left(1913\right)\left(653.5\right)-0.010672\left(1913\right)\\ =-457.9311\end{array}$

$\begin{array}{c}{K}_{2E_h}=0.00035\left(2185+\frac{-544}{2}\right)\left(681+\frac{-55}{2}\right)-0.00031\left(1741+\frac{502}{2}\right)\\ -0.010672\left(1741+\frac{502}{2}\right)\\ =0.00035\left(1913\right)\left(653.5\right)-0.00031\left(1992\right)-0.010672\left(1992\right)\\ =415.6748\end{array}$

$\begin{array}{c}{K}_{2I_h}=0.00031\left(1741+\frac{502}{2}\right)-0.071429\left(681+\frac{-55}{2}\right)-0.010672\left(681+\frac{-55}{2}\right)\\ =0.00031\left(1992\right)-0.071429\left(653.5\right)-0.010672\left(653.5\right)\\ =-53.0355\end{array}$

$\begin{array}{c}{K}_{2R_h}=0.071429\left(681+\frac{-55}{2}\right)-0.010672\left(1057+\frac{37}{2}\right)\\ =0.071429\left(653.5\right)-0.010672\left(1075.5\right)\\ =35.2011\end{array}$

$\begin{array}{c}{K}_{3S_h}=0.035341-0.00035\left(2185+\frac{-457.9311}{2}\right)\left(681+\frac{-53.0355}{2}\right)\\ -0.010672\left(2185+\frac{-457.9311}{2}\right)\\ =0.035341-0.00035\left(1956.03445\right)\left(654.48225\right)-0.010672\left(1956.03335\right)\\ =-468.9059\end{array}$

$\begin{array}{c}{K}_{3E_h}=0.00035\left(2185+\frac{-457.9311}{2}\right)\left(681+\frac{-53.0355}{2}\right)\\ -0.00031\left(1741+\frac{415.6748}{2}\right)-0.010672\left(174+\frac{415.6748}{2}\right)\\ =0.00035\left(1956.03445\right)\left(654.48225\right)-0.00031\left(1948.8374\right)\\ -0.010672\left(1948.8374\right)\\ =426.6643\end{array}$

$\begin{array}{c}{K}_{3I_h}=0.00031\left(1741+\frac{415.6748}{2}\right)-0.071429\left(681+\frac{-53.0355}{2}\right)\\ -0.010672\left(681+\frac{-53.0355}{2}\right)\\ =0.00031\left(1948.8374\right)-0.071429\left(654.48225\right)-0.010672\left(654.48225\right)\\ =-53.1295\end{array}$

$\begin{array}{c}{K}_{3R_h}=0.071429\left(681+\frac{-53.0355}{2}\right)-0.010672\left(1057+\frac{35.2011}{2}\right)\\ =0.071429\left(654.48225\right)-0.010672\left(1074.60055\right)\\ =35.2809\end{array}$

$\begin{array}{c}{K}_{4S_h}=0.035341-0.00035\left(2185-468.9059\right)\left(681-53.1295\right)\\ -0.010672\left(2185-468.9059\right)\\ =0.035341-0.00035\left(1716.0941\right)\left(627.8705\right)-0.010672\left(1716.0941\right)\\ =-395.3984\end{array}$

$\begin{array}{c}{K}_{4E_h}=0.00035\left(2185-468.9059\right)\left(681-53.1295\right)\\ -0.00031\left(1741+426.6643\right)-0.010672\left(1741+426.6643\right)\\ =0.00035\left(1716.0941\right)\left(627.8705\right)-0.00031\left(2167.6643\right)\\ -0.010672\left(2167.6643\right)\\ =353.3144\end{array}$

$\begin{array}{c}{K}_{4I_h}=0.00031\left(1741+426.6643\right)-0.071429\left(681-53.1295\right)\\ -0.010672\left(681-53.1295\right)\\ =0.00031\left(2167.6643\right)-0.071429\left(627.8705\right)-0.010672\left(627.8705\right)\\ =-50.8768\end{array}$

$\begin{array}{c}{K}_{4R_h}=0.071429\left(681-53.1295\right)-0.010672\left(1057+35.2809\right)\\ =0.071429\left(627.8705\right)-0.010672\left(1092.2809\right)\\ =33.1913\end{array}$

Hence, after a small change in time (or step in time) subsequent to the initial malaria incidence, we can assume new values of the model parameters, i.e.

$S_h\left(1\right)=2185+\left(-544\right)=1641\text{persons}$

$E_h\left(1\right)=1741+502=2243\text{persons}$

$I_h\left(1\right)=681+\left(-55\right)=626\text{persons}$

$R_h\left(1\right)=1057+37=1094\text{persons}$

From the above, we can see that the number of people exposed to the disease increases, which increases the likelihood of spread of the infection. Whereas the number of susceptible people and infected people decreases in number. We now have 1641 susceptible people, 2243 being exposed to the infection, 626 infected people and 1094 recovered people. These will now represent our new values that are used to continue modelling the spread of this particular disease over a period, recursively.

However, Table 2, below shows the calculated values of $S_h\left(t\right)$ , $E_h\left(t\right)$ , $I_h\left(t\right)$ and $R_h\left(t\right)$ modelled at time, t = 0, 1, 2, 3, and 4.

Table 2. Calculated values of SEIR model compartments based on recursion relation.

Time ($t_i$ )	Susceptible ($S_h$ )	Exposed ($E_h$ )	Infective ($I_h$ )	Recovered ($R_h$ )	Population Size ($N_h$ )	Effective Reproduction Number ($R_E$ )	Herd Immunity (v)
0	2185	1741	681	1057	5664	10.7065	0.9066
1	1641	2243	626	1094	5604	8.0409	0.8756
2	1183	2659	573	1129	5544	5.7967	0.8275
3	714	3086	520	1164	5484	3.4986	0.7142
4	318	3439	469	1197	5423	1.5582	0.3582

Table 2 shows that individuals susceptible to the disease decreases over time as the number of people exposed to the disease increases. However, there is a gradual decline in the number of infectious individuals perhaps due to the nature of the disease with an increase in the number of recovered persons. Hence, at t = 0,

$R_E=0.0049\times 2185=10.7065$ .

This is done for other values of $R_E$ at each time. Thus, $R_E$ of about 11 implies that one infected person can produce 11 new secondary cases of malaria incidence. This would spread faster in a population of susceptible individuals who are not yet vaccinated. However, if vaccination is done from the onset of disease incidence, then our $R_E$ reduces to about 2 at time t = 4, where an infected person will be able to produce 2 new secondary cases, perhaps as a result of reduction in susceptible individuals or an increase in vaccinated individuals.

Nevertheless, the dynamic trends for each of the four human compartments is seen from Figures 4(a)-(d) respectively.

From Figures 4(a)-(d), we see that as the proportion of susceptible people decreases, the proportion of people exposed to the disease increases, as shown respectively by blue and amber lines. However, there is a slow decline in the number of infectious individuals perhaps due to the nature and/or factors affecting disease spread and a slow increase in the number of recovered persons as seen from the red and green lines respectively. This means that the infection is being curtained in the population.

(a)

(b)

(c)

(d)

Figure 4. (a) Susceptible human compartment plot for Malaria over time; (b) Exposed human compartment plot for Malaria over time; (c) Infected human compartment plot for Malaria over time; (d) Recovered human compartment plot for Malaria over time.

Consequently, at the onset of disease incidence, our reproduction number for malaria outbreak from Equation (37), is given as,

$R=\frac{0.00035\times 1}{0.071429\times 1}=0.0049$

Hence, since our R value is 0.0049 < 1, this indicates that an epidemic is less likely to occur. Hence, the disease-free equilibrium for malaria transmission can be achieved as the number of infected individuals eventually declines in the population. Nevertheless, table 2 also shows that at the beginning of the Malaria incidence, our herd immunity threshold is 0.9066, this means that the government have to vaccinate about 91% of the population in order for the infection to be contained or to eradicate the disease from the population. However, after some change in time, our herd immunity threshold decreases to about 0.3582 at t = 4 which means that the number of infected individuals decreases as a result of increase in vaccination of the population at the onset of disease incidence. Thus, vaccination as well as good sanitary practice can help curb the spread of malaria infection among the populace in other to achieve herd immunity.

4. Conclusions

From the study, we had shown how compartmental modeling of the incidence of infectious diseases can be designed in compartments as seen from the SEIR-SEI model in Figure 2, the reproduction number determined and related to vaccination among study population.

Furthermore, while infectious diseases progress within populations both due to the behavior of the infectious agent and the population itself, Compartmental modelling on how they progress in an epidemic can inform how effective intervention will be (for example, social distancing or mass vaccination). There will always be new infectious viruses, bacteria, fungal and parasites with their corresponding outbreaks. However, with the aid of mathematical modelling such as compartmental SEIR-SEI model, epidemiologists, medical professionals and scientists will be able to develop a deeper understanding of the potentially harmful consequences of future diseases and to limit the effects they will have on the larger public.

Acknowledgements

I wish to express my profound gratitude to my research supervisors, Dr. U. P. Ogoke and Prof. E. C. Nduka for their guidance in the course of this research. In addition, my special thanks goes to my parents Obinna and Ngozi Enwere for their love, and to all my siblings for their support and encouragement.

Source of Funding

This research does not benefit from grant from any non-profit, public or commercial funding agency.

Consent for Publication

All authors contributed to the manuscript and consented to the publication of this research work.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	World Health Organization (2024) World Malaria Report 2023—Fact Sheets Overview.
[2]	Browne, C., Gulbudak, H. and Webb, G. (2015) Modeling Contact Tracing in Outbreaks with Application to Ebola. Journal of Theoretical Biology, 384, 33-49.[CrossRef] [PubMed]
[3]	Li, M.Y. and Muldowney, J.S. (1995) Global Stability for the SEIR Model in Epidemiology. Mathematical Biosciences, 125, 155-164.[CrossRef] [PubMed]
[4]	Liu, W., Hethcote, H.W. and Levin, S.A. (1987) Dynamical Behavior of Epidemiological Models with Nonlinear Incidence Rates. Journal of Mathematical Biology, 25, 359-380.[CrossRef] [PubMed]
[5]	Li, M.Y., Graef, J.R., Wang, L. and Karsai, J. (1999) Global Dynamics of a SEIR Model with Varying Total Population Size. Mathematical Biosciences, 160, 191-213.[CrossRef] [PubMed]
[6]	Jia, Z., Tang, G., Jin, Z., Dye, C., Vlas, S.J., Li, X., et al. (2008) Modeling the Impact of Immigration on the Epidemiology of Tuberculosis. Theoretical Population Biology, 73, 437-448.[CrossRef] [PubMed]
[7]	Bhunu, C.P., Garira, W., Mukandavire, Z. and Zimba, M. (2008) Tuberculosis Transmission Model with Chemoprophylaxis and Treatment. Bulletin of Mathematical Biology, 70, 1163-1191.[CrossRef] [PubMed]
[8]	Bishai, D., Johns, B., Nair, D., Nabyonga-Orem, J., Fiona-Makmot, B., Simons, E., et al. (2011) The Cost-Effectiveness of Supplementary Immunization Activities for Measles: A Stochastic Model for Uganda. The Journal of Infectious Diseases, 204, S107-S115.[CrossRef] [PubMed]
[9]	Wang, R., Jin, Z., Liu, Q., van de Koppel, J. and Alonso, D. (2012) A Simple Stochastic Model with Environmental Transmission Explains Multi-Year Periodicity in Outbreaks of Avian Flu. PLOS ONE, 7, e28873.[CrossRef] [PubMed]
[10]	Lekone, P.E. and Finkenstädt, B.F. (2006) Statistical Inference in a Stochastic Epidemic SEIR Model with Control Intervention: Ebola as a Case Study. Biometrics, 62, 1170-1177.[CrossRef] [PubMed]
[11]	Morton, A. and Finkenstädt, B.F. (2005) Discrete Time Modelling of Disease Incidence Time Series by Using Markov Chain Monte Carlo Methods. Journal of the Royal Statistical Society Series C: Applied Statistics, 54, 575-594.[CrossRef]
[12]	Macro Trends (2024) Nigeria Birth/Death Rate 1950-2024. https://www.macrotrends.net
[13]	Aakash, M., Gunasundari, C. and Qasem, M. (2023) Mathematical Modeling and Simulation of SEIR Model for COVID-19 Outbreak: A Case Study of Trivandrum. Frontiers in Applied Mathematics and Statistics, 9, Article ID: 1124897.[CrossRef]
[14]	Abigail, B. (2025) Numerical Analysis of the SEIR Model. University of Mary Washington.
[15]	Roberta, U. and Sudi, M. (2021) The Euler, Heun, and Fourth Order Runge-Kutta Solutions to SEIR Model for the Spread of Meningitis Disease. Mathline: Jurnal Matematika dan Pendidikan Matematika, 6, 140-153.[CrossRef]