Prerpints.org logo
Preprint
Article

A Mathematical Model for a Disease Outbreak Considering Waning-Immunity Class with Nonlinear Incidence and Recovery Rates

Altmetrics

Downloads

131

Views

34

Comments

0

A peer-reviewed article of this preprint also exists.

Submitted:

30 August 2023

Posted:

30 August 2023

You are already at the latest version

Alerts
Abstract
In this paper, we constructed and analyzed a modified SIR-type model for the epidemic problem which considering waning-immunity class in the population, with the nonlinear incidence and recovery rates. Depending on the basic reproduction ratio, we investigated conditions for both non-endemic and co-existing cases. The RouthHurwitz criteria is used to verify the local stability of equilibria, whereas for the global stability, the suitable Lyapunov function is selected to analyze the behavioral stability for each equilibria. Moreover, we studied the optimal control problem for this case. Numerically, we give some simulations to support our analitycal findings.
Keywords: 
Subject: Public Health and Healthcare  -   Health Policy and Services

1. Introduction

The behavior of the spread of disease can be studied mathematically through a compartment model. One model that is often used is the SIR (Susceptible-Infectious-Recovered) model. This model extends the SI model, which previously did not consider the recovery class. In its development, the SIR model has been modified by considering various aspects. Next, various factors that influence disease-spreading behavior have been reasons for some researchers to consider them in their work. For examples, quarantine [3], treatment [4], diffusion [5,6], vaccination [8], delayed [9], combination from delayed, vaccination and treatment [7], social distancing [39,40], relapse and media impact [10], and others [37,38,41,42,43].
The various studies mentioned above still apply a bilinear infection rate. This differs from some other work on epidemic modeling, which has applied a nonlinear incidence rate form. Regarding the application of the nonlinear incidence rate within the SIR model, several work has been done previously. For examples, Jin et al. in [33], Zhang et al. in [24], McCluskey in [25], Enatsu and Nakata in [26], Li et al. in [34], Chen and Zhao in [27], Ammi et al. in [28], Koufi et al. in [31], and others (see [29,30,32]).
Hereinafter obviousness that the treatment is a main strategic to handle diseases outbreak. In the mathematical epidemiology, generally it was assumed to be the recovery rate in consequence of treatment and expressed in the form: T ( I ) = r I , I ≤ 0 . Nevertheless, in its application the formula is ineffective if the availability of healthcare resources such as vaccines, drugs, hospital beds, etc., are insufficient amount. However, the supply of those patient amenities in hospital care is always limited contextually. Further, for this reason, several researchers have applied nonlinear forms of treatment to their work. For examples, Zhang and Liu [13] used nonlinear treatment form is r I 1 + α I , where r > 0 is represent signifies the cure rate and α ≥ 0 is indicates the significant impacts of the infection on the people who are being delayed for treatment. Next, Rajasekar et al. in [15] has analyzed the model of Zhang and Liu in [13] using a stochastic approach. Ghosh et al. in [14] proposed and analyzed a model with inhibitory effect and saturated treatment function, where nonlinear treatment form is r u 2 I 1 + b u 2 I ( u 2 represent optimal treatment strategies to minimize the infected cases).
Due to the health planning to the public service, a standard to estimate resource availibility was fixed by the World Health Organization (WHO) Statistical Information System i.e., available hospital beds per 10,000 population [1,11,12,36]. By considering this situation, next Shan and Zhu [36], Li and Zhang [12], Cui et el. [11], as well as Alshammari and Khan [1] further modified the proportinal of treatment rate, which the capacity of hospital was considered, where T ( I ) = α 0 + α 1 − α 0 m m + I . Parameter α 0 and α 1 are, respectively, the minimum and maximum per capita recovery rates. A constant m is represents the measure of the quantity of available hospital beds.
The strategy combinations between the government action and hospitalization situation need to be incorporated to handle a better of disease outbreak. Furthermore, the combination of nonlinear incidence and treatment rates has become interesting and has been carried out by several researchers [11,12,14,36]. Recently, Alshammari and Khan [1] set incidence function in the form of β S I k + I and treatment function in the form of α 0 + α 1 − α 0 m m + I . They found that for limited hospital beds led to the system undergoes bifurcation. In addition, by increasing the value of the intervention k has a significant impact on reducing the case quantity in the infection class.
On the other side, the behaviors of immunity condition on the host Individuals plays a vital role in the transmission of a diseases. Individuals which have just healed from disease get any immunity, however, its standard of effectiveness is not necessarily constant respect to time. An individual’s immunity may decay as time goes by, it could be the individual will be reinfection due to exposure to infectious agent [23]. For examples: measles [19,21,22], mumps [20,21,22], rubella [21,22] and COVID-19 [2] are of this type. Thus, waning immunity in individuals ought to serious attention because it has the opportunity to increase the number of cases of reinfection.
Motivated by these various studies, we propose an epidemic model by considering waning immunity with nonlinear incidence and recovery rates. This model extended from a model by Alshammari and Khan in [1]. The paper is organized as follows. In Section 2 we establish mathematical model for the current epidemic problem. We analyze the properties of the model in Section 3. These analysis include: positivity and invariant region, the existences of the non-endemic and co-existing equilibria, as well as determine a basic reproduction ratio. Furthermore, we find conditions for the local and global behaviours of the each equilibria. Next, The optimal control problem is shown in Section 4. Section 5, we conduct numerical simulations to confirm our analytical findings. These simulations are performed for the various of the parameters values. Finally, our work deduced in Section 6.

2. Model Formulation

Thoroughly, the population is M ( t ) , and next, it will be classified into four categories: susceptible S ( t ) , infectious I ( t ) , recovered R ( t ) , and susceptible that previously infected W ( t ) humans, where M = S + I + R + W (as shown in Figure 1). Individuals move from compartment S to I as well as W to I as with transmission rates: β 1 I S k + I and β 2 I W k + I . We have four sub–populations, therefore a nonlinear dynamical system consisting of four nonlinear differential equations established is as follows:
d S d t = b − β 1 I S k + I − μ S d I d t = β 1 I S k + I + β 2 I W k + I − α 0 + α 1 − α 0 m m + I I − γ + μ I d R d t = α 0 + α 1 − α 0 m m + I I − μ R − ξ R d W d t = ξ R − β 2 I W k + I − μ W
Table 1. Parameters Description
Table 1. Parameters Description
Parameters Descriptions Units
b The birth rate P e o p l e × T i m e − 1
β 1 The probability of transmission from S to I ( P e o p l e × T i m e ) − 1
β 2 The probability of transmission from E to I ( P e o p l e × T i m e ) − 1
μ The natural death rate T i m e − 1
γ The death due to disease T i m e − 1
ξ The probability rate of waning immune of people T i m e − 1
α 0 The minimum per capita recovery rate T i m e − 1
α 1 The maximum per capita recovery rate T i m e − 1
k The intervention levels P e o p l e
m The impact of the number of hospital beds on the ourbreak P e o p l e

3. Mathematical Analysis

3.1. Positivity of Solutions

Lemma 1. 
With the non-negative initial values, say χ ( 0 ) and χ ( 0 ) = ( S , I , R , W ) representing the state variables, each of the solutions curves of the model (1) will be non-negative for all t > 0 . Additionally,
lim x → ∞ sup M ( t ) ≤ b μ .
Proof. 
Let t 1 = sup t > 0 : χ ( t ) > 0 , so t 1 > 0 , and it follows by the first equation of system (1) that:
d S ( t ) d t = b − β 1 I ( t ) S ( t ) k + I ( t ) − μ S ( t ) .
It can be re-written as:
d S ( t ) d t exp ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ + μ t + S ( t ) β 1 I ( τ ) k + I ( τ ) + μ exp ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ + μ t = b exp ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ + μ t .
Therefore,
d d t S ( t ) exp ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ + μ t = b exp ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ + μ t .
Hence,
S ( t ) exp ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ + μ t − S ( 0 ) = b exp ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ + μ t .
So,
S ( t ) = S ( 0 ) exp − ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ − μ t + exp − ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ − μ t + b exp ∫ 0 t β 1 I ( τ ) k + I ( τ ) d τ + μ t > 0 .
For the rest of the equations, we can consistently determine the positivity via this method. So, we are able to say that the solution of the model given by (1) is non-negative for every time t > 0 . For the rest of the part of the proof in the considered Theopositivity, we have 0 < χ ( 0 ) ≤ N ( t ) , and summing the equations involve in the system (1), we arrive at the following expression,
d M d t = b − μ M − γ I ≤ b − μ M , thus lim x → ∞ sup M ( t ) ≤ b μ .
 □

3.2. Invariant Region

For the biological significance of our model (1), we will show that the variables and the parameters are non-negative for all time t ≥ 0 and analyse the model (1) in a suitable feasible region Φ .
Lemma 2. 
The feasible region Φ defined by:
Φ = ( S , I , R , W ) ∈ R + 4 : S + I + R + W ≤ b μ ,
with initial condition S ( 0 ) ≥ 0 , I ( 0 ) ≥ 0 , R ( 0 ) ≥ 0 , W ( 0 ) ≥ 0 is positively invariant for system (1).
Proof. 
Adding the equations of system (1), we have:
d M d t ≤ b − μ M .
It follows that 0 ≤ M ( t ) ≤ b μ + M ( 0 ) e − μ t , where M ( 0 ) denotes the initial values of the total population. Thus 0 ≤ M ( t ) ≤ b μ , as t → ∞ . So, the region:
Φ = ( S , I , R , W ) ∈ R + 4 : S + I + R + W ≤ b μ ,
is a positively invariant set for system (1). □

3.3. Non-endemic Equilibria and Basic Reproduction Ratio

The non-endemic equilibria of the model (1) is obtained by setting I = 0 , and substituting it into (1) to obtain:
E 0 = S 0 , I 0 , R 0 , W 0 = b μ , 0 , 0 , 0 , 0 .

3.4. The Basic Reproduction Ratio

Here, we have the following matrix of new infection F , and the matrix of transfer V . We can write
F = β 1 I S k + I + β 2 I W k + I + α 0 m I m + I 0 0 0 ,
V = α 0 + α 1 m m + I I + γ + μ I − b + β 1 I S k + I + μ S − ξ R + β 2 I W k + I + μ W − α 0 + α 1 − α 0 m m + I I + μ R + ξ R .
The Jacobians of the above matrices at the non-endemic equilibria P 0 respectively are respectively given by
F = α 0 + b β 1 k μ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ,
V = α 0 + α 1 + γ + μ 0 0 0 b β 1 k μ μ 0 0 0 0 μ 0 − α 1 0 0 μ + ξ .
Next, we get the eigenvalues of matrix F V − 1 are λ 1 , 2 , 3 = 0 and λ 4 = α 0 k μ + b β 1 k μ ( α 0 + α 1 + γ + μ ) . Moreover, the reproduction ratio is
R 0 = ρ ( F V − 1 ) = α 0 k μ + b β 1 k μ ( α 0 + α 1 + γ + μ ) .

3.5. The Co-existing Equilibria

System (1) has a unique co-existing equilibria in the interior of Φ that is given by
E * = S * , I * , R * , W *
where
S * = b ( k + I * ) k μ + ( β 1 + μ ) I * , R * = α 0 I * ( m + I * ) + m ( α 1 − α 0 ) I * ( ξ + μ ) ( m + I * ) , W * = ξ ( k + I * ) R * k μ + ( β 2 + μ ) I *
Next, we substitution S * , R * and W * into the second equation in (1) so yield a cubic equation in I * as follow,
C 1 ( I * ) 3 + C 2 ( I * ) 2 + C 3 ( I * ) + C 4 = 0 ,
where
C 1 = ( β 1 + μ ) ( β 2 + μ ) ( ξ + μ ) α 1 R 0 α 0 k μ + b β 1 α 1 k μ − R 0 C 2 = ( ξ + μ ) ( β 1 + μ ) ( β 2 + μ ) α 0 R 0 α 0 k μ + b β 1 α 0 k μ − R 0 m + ( ξ + μ ) α 1 R 0 α 0 k μ + b β 1 α 1 k μ − R 0 k μ [ ( β 1 + μ ) + ( β 2 + μ ) ] , C 3 = ( ξ + μ ) k μ α 1 R 0 α 0 k μ + b β 1 α 1 k μ − R 0 k μ + ( ξ + μ ) k μ α 0 R 0 α 0 k μ + b β 1 α 0 k μ − R 0 m [ ( β 1 + μ ) + ( β 2 + μ ) ] − [ ( β 1 + μ ) ( ξ m α 1 β 2 ) + b β 1 ( β 2 + μ ) ( ξ + μ ) ] , C 4 = ( ξ + μ ) ( α 0 k μ + b β 1 ) m α 0 R 0 1 α 0 − 1 α 0 + α 1 + γ + μ − k μ b β 1 ( ξ + μ ) + ξ m α 1 β 2 .
which for the polynomial (5) have the roots I = 0 or I 1 , 2 = I 1 , 2 * for C 4 = 0 , which can be written by
I 1 , 2 * = − C 2 ± C 2 2 − 4 C 1 C 3 2 C 1 .
From (5), it can be seen that C j < 0 , j = 1 , 2 , 3 ⇔ R 0 > 1 Hence, we have the following theorem regarding the existence of the co-existing equilibria when R 0 > 1 .

3.6. Local Stability of Non-Endemic Equilibria

To investigate the local stability of all equilibria, we linearize (1), and it yield a Jacobian matrix,
J = ( − β 1 I k + I + μ − k β 1 S ( k + I ) 2 0 0 β 1 I k + I k β 1 S ( k + I ) 2 + k β 2 W ( k + I ) 2 − ( α 0 + γ + μ ) − ( α 0 − α 1 ) m 2 ( m + I ) 2 0 β 2 I k + I 0 α 0 + ( α 0 − α 1 ) m 2 ( m + I ) 2 − ( ξ + μ ) 0 0 − k β 2 S ( k + I ) 2 ξ − β 2 I k + I + μ )

3.7. Local stability of non-endemic equilibria

Theorem 1. 
The non-endemic equilibria E 0 is locally asymptotically stable for R 0 < 1 and unstable for R 0 > 1 .
Proof. 
From Jacobian matrix (7) at point E 0 , it is yields eigenvalues λ 1 , 2 = − μ . In addition, we get a characteristic equation,
f 01 λ 2 + f 11 λ + f 21 = 0 ,
where
f 01 = 1 , f 11 = 1 R 0 ( ξ + μ ) R 0 + α 1 α 0 k μ + b β 1 α 1 k μ − R 0 , f 21 = ( ξ + μ ) α 1 R 0 α 0 k μ + b β 1 α 1 k μ − R 0 .
Next, choose α 0 k μ + b β 1 = α 1 k μ . Moreover, if R 0 < 1 then the values of f 11 and f 21 are positive. Thus, E 0 is locally asymptotically stable whenever conditions f 01 > 0 , f 11 > 0 , f 21 > 0 , and f 11 f 21 > 0 is fullfiled.  □

3.8. Local stability of co-existing equilibria

Theorem 2. 
The co-existing equilibria E * is locally asymptotically stable whenever it exists.
Proof. 
From Jacobian matrix (7) at point E * , we get a characteristic equation,
f 02 λ 4 + f 12 λ 3 + f 22 λ 2 + f 32 λ + f 42 = 0 ,
where
f 02 = 1 , f 12 = Ψ 1 + ξ + μ − ( Ψ 2 + Ψ 3 ) , f 22 = Θ 1 Θ 2 + Θ 4 Θ 5 + ( ξ + μ ) ( Ψ 1 − ( Ψ 2 + Ψ 3 ) ) + Ψ 2 Ψ 3 − Ψ 1 ( Ψ 2 + Ψ 3 ) , f 32 = ( ξ + μ ) Ψ 2 Ψ 3 + Θ 1 Θ 2 + Θ 4 Θ 5 − Ψ 1 ( Ψ 2 + Ψ 3 ) + Ψ 1 Θ 4 Θ 5 + Ψ 3 ( Ψ 1 Ψ 2 − Θ 1 Θ 2 ) − ξ Θ 5 ( Θ 3 + α 0 ) , f 42 = ( ξ + μ ) ( ( Ψ 1 Ψ 2 − Θ 1 Θ 2 ) Ψ 3 + Ψ 1 Θ 4 Θ 5 ) − ξ Ψ 1 Θ 5 ( Θ 3 + α 0 ) , Ψ 1 = Θ 1 + μ , Ψ 2 = Θ 2 + Θ 4 − ( Θ 3 + α 0 + γ + μ ) , Ψ 3 = Θ 5 + μ , Θ 1 = β 1 I * k + I * , Θ 2 = k β 1 S * ( k + I * ) 2 , Θ 3 = ( α 0 − α 1 ) m 2 ( m + I * ) 2 , Θ 4 = k β 2 W * ( k + I * ) 2 , Θ 5 = β 2 I * k + I * .
From (9), we obtain that λ j , j = 1 , 2 , 3 , 4 will be negative if f i 2 > 0 , i = 1 , 2 , 3 , 4 , R 0 > 1 , f 12 f 22 > f 32 f 02 , and f 12 f 22 f 32 > f 12 2 f 42 + f 02 f 32 2 . If this condition holds, then based on the Routh-Hurwitz criterion, point E * is locally asymptotically stable.  □
From Theorems 1 and 2, we note that if the value of the basic reproduction ratio is less than one, then the non-endemic equilibria E 0 will be asymptotically stable. Otherwise, the co-existing equilibria E * exist and are stable whenever the reproduction ratio is greater than one. Furthermore, we remark that when R 0 < 1 the co-existing equilibria E * does not exist. In this case, R 0 is a threshold. On the contrary, if R 0 > 1 the E 0 becomes unstable and E * exist. If the Routh-Hurwitz criterion of equation (9) holds, then E * is asymptotically stable.
Corollary 1. 
If R 0 > 1 , f i 2 > 0 , i = 1 , 2 , 3 , 4 , f 12 f 22 > f 32 f 02 , and f 12 f 22 f 32 > f 12 2 f 42 + f 02 f 32 2 then the system (1) undergoes a forward bifurcation at R 0 = 1 from E 0 to E * .

3.9. Global stability of non-endemic equilibria

Theorem 3. 
The non-endemic equilibria E 0 is globally asymptotically stable for R 0 < 1 and unstable for R 0 > 1 .
Proof. 
Refer to global proving by the various works in [16,17,18,35], define a Lyapunov function
L 0 = S − S 0 − S 0 ln S S 0 + I + R + W
Differentiating with respect to time yields
d L 0 d t = 1 − S 0 S d S d t + d I d t + d R d t + d W d t
Next, from (1) we get
d L 0 d t = 1 − S 0 S b − β 1 I S k + I − μ S + β 1 I S k + I + β 2 I W k + I − α 0 + α 1 − α 0 m m + I I − γ + μ I + α 0 + α 1 − α 0 m m + I I − ( μ + ξ ) R + ξ R − β 2 I W k + I − μ W = μ S 1 0 2 − S S 0 − S 0 S + ( b 1 S 0 − k ( γ + μ ) ) I − ( γ + μ ) I 2 − μ ( R + W ) ≤ μ S 0 2 − S S 0 − S 0 S + R 0 − 1 + α 0 α 0 + α 1 + γ + μ I − ( γ + μ ) I 2 − μ ( R + W ) ≤ μ S 1 0 2 − S S 0 − S 0 S + ( R 0 − 1 ) I − ( γ + μ ) I 2 − μ ( R + W )
The value of d L 0 d t < 0 if R 0 < 1 . Using the relation between the geometric means and arithmetic means, we confirm that d L 0 d t ≤ 0 and fulfilled the equality only at E 0 . Thus, the non-endemic equilibria E 0 is globally asymptotically stable if R 0 < 1 .  □

3.10. Global stability of co-existing equilibria

Theorem 4. 
The co-existing equilibria E * is globally asymptotically stable whenever it exists.
Proof. 
Refer to global proving by the various works in [16,17,18,35], define a Lyapunov function
L * = S − S * − S * ln S S * + I − I * − I * ln I I * + R − R * − R * ln R R * + W − W * − W * ln W W * .
Differentiating with respect to time yields
d L * d t = 1 − S * S d S d t + 1 − I * I d I d t + 1 − R * R d R d t + 1 − W * W d W d t .
Next, from (1) we get
d L * d t = 1 − S * S b − β 1 I S k + I − μ S + 1 − I * I β 1 I S k + I + β 2 I W k + I − α 0 + α 1 − α 0 m m + I I − γ + μ I + 1 − R * R α 0 + α 1 − α 0 m m + I I − μ R − ξ R + 1 − W * W ξ R − β 2 I W k + I − μ W = μ S * 2 − S S * − S * S + β 1 I * S * k + I * 2 − I I * − S S * k + I * k + I − S * S + β 2 I * W * k + I * 2 − I I * − W W * k + I * k + I + C 1 − C 2
where
C 1 = β 1 I S * k + I + β 2 I * W k + I * + I W * k + I + ξ R + R * W + α 0 + α 1 − α 0 m m + I I + α 0 + α 1 − α 0 m m + I * I * , C 2 = 2 β 2 I * W * k + I * + ξ R * W * + W * R W + α 0 + α 1 − α 0 m m + I * I * R R * + α 0 + α 1 − α 0 m m + I I R * R .
Next, since a relation between geometric means and arithmetic means is applied, we claim that d L * d t ≤ 0 . This condition holds only at E * . In consequence, a co-existing point E * is globally asymptotically stable.  □

4. Optimal Control Problem

To control the spread of the disease, we set to suppress the number of transitioned individuals from recovered into susceptible again. Procedures to suppress that transition phenomenon can be conducted through educational efforts, including outpatient; consultancy; or health campaigns. To reduce the number of the W population with the minimum cost, we constructed the dynamical model in system (1) by adding parameters for the control of educational efforts, namely u ( t ) . Then, we obtain:
d S ( t ) d t = b − β 1 I ( t ) S ( t ) k + I ( t ) − μ S ( t ) d I ( t ) d t = β 1 S ( t ) + β 2 W ( t ) k + I ( t ) I ( t ) − α 0 + α 1 − α 0 m m + I ( t ) I ( t ) − γ + μ I ( t ) d R ( t ) d t = α 0 + α 1 − α 0 m m + I ( t ) I ( t ) − μ R ( t ) − ( 1 − u ( t ) ) ξ R ( t ) d W ( t ) d t = ( 1 − u ( t ) ) ξ R ( t ) − β 2 I ( t ) W ( t ) k + I ( t ) − μ W ( t )
We determine educational effort as an intervention with the objective functions as follows:
J ( u ) = min ∫ 0 T [ C 1 W ( t ) + C 2 u 2 ( t ) ] d t
where T is the final time, and C 1 , C 2 is the weight constant and cost of the conducted interventions. The control problem formed from the constraint equation as the equation of each compartment (10) and the objective function (11) can form a Hamiltonian equation related to the optimal control problem. Therefore, we use Pontryagin’s maximum principle for the optimal control ( u * ) ∈ U satisfying (10), such that the associated pseudo-Hamiltonian is
H = C 1 W ( t ) + C 2 u 2 ( t ) + L 1 b − β 1 I ( t ) S ( t ) k + I ( t ) − μ S ( t ) + L 2 β 1 S ( t ) + β 2 W ( t ) k + I ( t ) I ( t ) − α 0 + α 1 − α 0 m m + I ( t ) I ( t ) − γ + μ I ( t ) + L 3 α 0 + α 1 − α 0 m m + I ( t ) I ( t ) − μ R ( t ) − ( 1 − u ( t ) ) ξ R ( t ) + L 4 ( 1 − u ( t ) ) ξ R ( t ) − β 2 I ( t ) W ( t ) k + I ( t ) − μ W ( t )
where L i , i = 1 , . . . , 4 are adjoint variables satisfied i.e.
L ˙ 1 = − L 1 − β 1 I ( t ) k + I ( t ) − μ − L 2 β 1 I ( t ) k + I ( t ) L ˙ 2 = − L 1 − β 1 S ( t ) k + I ( t ) + β 1 I ( t ) S ( t ) k + I ( t ) 2 − L 2 β 1 S ( t ) k + I ( t ) − β 1 I ( t ) S ( t ) k + I ( t ) 2 + β 2 W ( t ) k + I ( t ) − L 2 β 2 I ( t ) W ( t ) k + I ( t ) 2 + ( α 1 − α 0 ) m I ( t ) m + I ( t ) 2 − α 0 − ( α 1 − α 0 ) m m + I ( t ) − γ − μ − L 3 − ( α 1 − α 0 ) m I ( t ) m + I ( t ) 2 + α 0 + ( α 1 − α 0 ) m m + I ( t ) − L 4 − β 2 W ( t ) k + I ( t ) + β 2 I ( t ) W ( t ) k + I ( t ) 2 L ˙ 3 = − L 3 − μ − 1 − u ( t ) ξ − L 4 1 − u ( t ) ξ L ˙ 4 = − C 1 − L 2 β 2 I ( t ) k + I ( t ) − L 4 − β 2 I ( t ) k + I ( t ) − μ
where the final condition L i ( T ) = 0 for i = 1 , . . . , 4 . The necessary and sufficient optimal condition satisfying are obtained, which in turn gives the optimal control
u * = min 1 , max 0 , 1 2 ξ R ( t ) L 3 − L 4 C 2 .

5. Numerical Simulation

This numerical simulation is designed to support the results of the analysis discussed in the previous section. The dynamical population will be compared when some parameter values change. We applied Runge–Kutta Fourth–Order to solve model (1) using the parameter values in Table 2 and initial values in Table 3.
Figure 2 shown that E * stable if β 1 = 0.01 . The population in each compartment converges to a certain amount. This indicates that individuals still exist in the population at a certain time interval. Clearly, it is different if β 1 = 0.0018 , which E 0 stable. This case is depicted in Figure 3. The population in each compartment I, R, and W converge to zero. This situation showed that the disease is free in the system. The behavior of the system (1) related to the changes of parameter β 1 from 0.0018 to 0.01. If β 1 > β 1 * then E * stable, meanwhile E 0 unstable precisely when β 1 < β 1 * , where β 1 * = 0.001804347826 . This indicates that the system (1) undergoes a forward bifurcation. This case has been illustrated in the Figure 4.
Figure 5a shown that the increase in infection ( β 1 ) causes R 0 to also enlarge. On the other hand, increased intervention level k can lead to minimizing R 0 . Thereby, intervention is an important factor in reducing the number of infection cases. Next, Figure 5b shows that the smaller α 0 and α 1 causes R 0 to increase. This illustrates that a lower recovery rate significantly impacts the increase in the number of infections. The condition of waning immunity experienced by individuals who have recovered will certainly affect the dynamics of the system. An increase in the value of ξ causes the population in R to decrease. Meanwhile, the population in W increasing. This behavior is plotted in Figure 6.

5.1. Numerical Simulation on Optimal Control Problem

Respective to the objective function on Equation 11, we aim to minimize the number of infected, uneducated injector populations, treatment, and educational efforts. The optimal function graph is shown in the following figure.
Based on Figure 8, it is clear that the system with applied optimal control successfully suppresses the number W population since they were staying on the R population after healing from infection. The number of W population stays at the bottom along the control applied, but when at the end, the W population increase. Referring to Figure 7, it happens since the control value seems instantaneously to zero. Nevertheless, the optimal control success to reach the aim of the objective function is to minimize the W population and cost of educational effort.

6. Conclusions

We have formulated a SIRW–type model for an epidemic problem considering waning immunity class with nonlinear incidence and recovery rates. The model is class-based in the form of a differential equations system, where the population is divided into Susceptible ( S ) , Infectious ( I ) , Recovered ( R ) , Susceptible that previously infected ( W ) . Using a method, namely a next-generation matrix, we obtained the Basic Reproduction Ratio ( R 0 ) , which is a threshold to control the transmission of disease. Then the simulation results show that the waning immunity factor and half-saturated infection can affect disease transmission. The period of waning immunity and the grade of saturated incidence can spread of disease slowly. The results obtained can be used to reference early prevention of the spread of disease as long as it has similar spreading behavior.

Author Contributions

Conceptualization, N.A. and L.K.B.; methodology, N.A. and L.K.B.; software, N.A., L.K.B., F.I., and S.T.T; validation, N.A. and L.K.B.; formal analysis, N.A. and L.K.B.; investigation, N.A. and L.K.B.; resources, N.A., L.K.B., F.I., and S.T.T; writing—original draft preparation, N.A., L.K.B., F.I., and S.T.T; writing—review and editing, N.A., L.K.B., F.I., and S.T.T; visualization, N.A., L.K.B., F.I., and S.T.T; supervision, N.A.; project administration, N.A., L.K.B., F.I., and S.T.T; funding acquisition, N.A. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Universitas Padjadjaran, Indonesia, via Hibah Riset Data Pustaka dan Daring Universitas Padjadjaran, No. 1549/UN6.3.1/PT.00/2023.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data sharing not applicable.

Acknowledgments

The authors would like to thank Universitas Padjadjaran to support our project. Also, this project is part of the post doctoral research conducted by LKB.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Proof of the numerical analysis of the non-endemic equilibria

Here, numerically, the Routh-Hurwitz criterion is implemented to assert the Theorem 1 as follows ( β 1 = 0.0018 ):
f 01 = 1 > 0 , f 11 = 0.026 > 0 , f 21 = 0 : 000025 > 0 .
In addition to R 0 = 0.9983739837 < 1 .

Appendix B. Proof of the numerical analysis of the co-existing equilibria

The similar method in the Appendix A is also used to assert the Theorem 2 as follows ( β 1 = 0.01 ):
f 02 = 1 > 0 , f 12 = 0.3234529156 > 0 , f 22 = 0.01496557001 > 0 , f 32 = 0.0002090709850 > 0 , f 42 = 8.954418034 × 10 − 7 > 0 , f 12 f 22 − f 32 f 02 = 0.004631586268 > 0 , f 12 f 22 f 32 − f 12 2 f 42 + f 02 f 32 2 = 8.746475801 × 10 − 7 > 0 .
In addition to R 0 = 4.065040650 > 1 .

References

  1. Alshammari, F.S.; Khan, M.A. Dynamic behaviors of a modified SIR model with nonlinear incidence and recovery rates. Alexandria Engineering Journal 2021, 60, 2997–3005. [Google Scholar] [CrossRef]
  2. Anggriani, N.; Ndii, M. Z.; Amelia, R.; Suryaningrat, W.; Pratama, M. A. A. A mathematical COVID-19 model considering asymptomatic and symptomatic classes with waning immunity. Alexandria Engineering Journal 2022, 61, 113–124. [Google Scholar] [CrossRef]
  3. Beay, L. K. Modelling the effects of treatment and quarantine on measles. AIP Conference Proceedings 2018, 1937, 020004. [Google Scholar]
  4. Rafiq, M.; Ali, J.; Riaz, M.B.; Awrejcewicz, J. Numerical analysis of a bi-modal covid-19 SITR model. Alexandria Engineering Journal 2022, 61, 227–235. [Google Scholar] [CrossRef]
  5. Gai, C.; Iron, D.; Kolokolnikov, T. Localized outbreaks in an S-I-R model with diffusion. Journal of Mathematical Biology 2020, 80, 1389–1411. [Google Scholar] [CrossRef]
  6. Singh, H.P.; Bhatia, S.K.; Jain,R.; Bahri, Y. A Study on the Effect of Optimal Control Strategies: An SIR Model with Delayed Logistic Growth, In: Sharma, T.K., Ahn, C.W., Verma, O.P., Panigrahi, B.K. (eds) Soft Computing: Theories and Applications. Advances in Intelligent Systems and Computing, Springer, Singapore, 2021; 1381.
  7. Elazzouzi, E.; Alaoui, A.L.; Tilioua, M.; Tridane, A. Global stability analysis for a generalized delayed SIR model with vaccination and treatment. Advances in Difference Equations 2019, 532, 1–19. [Google Scholar] [CrossRef] [PubMed]
  8. d’Onofrio, A.; Manfredi, P. Bifurcation Thresholds in an SIR Model with Information-Dependent Vaccination. Mathematical Modelling of Natural Phenomena 2007, 2(1), 26–43. [Google Scholar] [CrossRef]
  9. Ebraheem, H.K.; Alkhateeb, N.; Badran, H.; Sultan, E. Delayed Dynamics of SIR Model for COVID-19. Open Journal of Modelling and Simulation 2021, 9, 146–158. [Google Scholar] [CrossRef]
  10. Rajasekar, S.P.; Quanxin Zhu, M. Higher order stochastically perturbed SIRS epidemic model with relapse and media impact. Open Journal of Modelling and Simulation 2022, 45(2), 843–863. [Google Scholar]
  11. Cui, Q.; Qui, Z.; Liu, W.; Hu, Z. Complex dynamics of an SIR epidemic model with nonlinear saturate incidence and recovery rate. Entropy 2017, 19(7), 305. [Google Scholar] [CrossRef]
  12. Li, G-H.; Zhang, Y.X. Dynamic behaviors of a modified SIR model in epidemic diseases using nonlinear incidence and recovery rates. Plos One 2017, 12(4), e0175789. [Google Scholar] [CrossRef] [PubMed]
  13. Zhang, X.; Liu, X. Backward bifurcation of an epidemic model with saturated treatment function. Journal of mathematical analysis and applications 2008, 348, 433–443. [Google Scholar] [CrossRef]
  14. Ghosh, J.K.; Ghosh, U.; Biswas, M.H.A.; Sarkar, S. Qualitative Analysis and Optimal Control Strategy of an SIR Model with Saturated Incidence and Treatment. Differential Equations and Dynamical Systems 2023, 31, 53–67. [Google Scholar] [CrossRef]
  15. Rajasekar, S.P.; Pitchaimani, M.; Zhu, Q. Dynamic threshold probe of stochastic SIR model with saturated incidence rate and saturated treatment function. Physica A 2019, 535, 122300. [Google Scholar] [CrossRef]
  16. Beay, L.K.; Anggriani, N. Dynamical analysis of a modified epidemic model with saturated incidence rate and incomplete treatment. Axioms 2022, 11(6), 256: 1–21. [Google Scholar] [CrossRef]
  17. Anggriani, N.; Beay, L.K. Modeling of COVID-19 spread with self-isolation at home and hospitalized classes. Results in Physics 2022, 36, 105378. [Google Scholar] [CrossRef] [PubMed]
  18. Samui, P.; Mondal, J.; Khajanchi, S. A mathematical model for COVID-19 transmission dynamics with a case study of India. Chaos, Solitons and Fractals 2020, 140, 110173. [Google Scholar] [CrossRef]
  19. Perry, R.T.; Halsey, N.A. The Clinical Significance of Measles: A Review. The Journal of Infectious Diseases 2004, 189 (Suppl. 1), S4–S16. [Google Scholar]
  20. Conly, J.M.; Johnston, B.L. Is mumps making a comeback? Canadian Journal of Infectious Disease and Medical Microbiology 2007, 18, 7–9. [Google Scholar] [CrossRef] [PubMed]
  21. Bankamp, B.; Hickman, C.; Icenogle, J.P.; Rota, P.A. Successes and challenges for preventing measles, mumps and rubella by vaccination. Current Opinion in Virology 2019, 34, 110–116. [Google Scholar] [CrossRef] [PubMed]
  22. Yang, L.; Grenfell, B.T.; Mina, M.J. Waning immunity and re-emergence of measles and mumps in the vaccine era. Current Opinion in Virology 2020, 40, 48–54. [Google Scholar] [CrossRef] [PubMed]
  23. Okuwa, K.; Inaba, H.; Kuniya, T. An age-structured epidemic model with boosting and waning of immune status. Mathematical Biosciences and Engineering 2021, 18(5), 5707–5736. [Google Scholar] [CrossRef] [PubMed]
  24. Zhang, J-Z.; Jin, Z.; Liu, Q-X.; Zhang, Z-Y. Analysis of a delayed SIR model with nonlinear incidence rate. Discrete Dynamics in Nature and Society 2008, 636153. [Google Scholar] [CrossRef]
  25. McCluskey, C.C. Global stability of an SIR epidemic model with delay and general nonlinear incidence. Mathematical Biosciences and Engineering 2010, 7(4), 837–850. [Google Scholar] [PubMed]
  26. Enatsu, Y.; Nakata, Y. Stability and bifurcation analysis of epidemic models with saturated incidence rates: An application to a nonmonotone incidence rate. Mathematical Biosciences and Engineering 2014, 11(4), 785–805. [Google Scholar] [CrossRef]
  27. Chen, Y.; Zhao, W. Dynamical analysis of a stochastic SIRS epidemic model with saturating contact rate. Mathematical Biosciences and Engineering 2020, 17(5), 5925–5943. [Google Scholar] [CrossRef] [PubMed]
  28. Ammi, M.R.S.; Tahiri, M.; Torres, D.F.M. Global stability of a caputo fractional SIRS model with general incidence rate. Mathematics in Computer Science 2021, 15, 91–105. [Google Scholar] [CrossRef]
  29. Thirthar, A.A.; Naji, R.K.; Bozkurt, F.; Yousef, A. Modeling and analysis of an SI1I2R epidemic model with nonlinear incidence and general recovery functions of I1. Chaos, Solitons and Fractals 2021, 145, 110746. [Google Scholar] [CrossRef]
  30. Sun, R. Global stability of the endemic equilibrium of multigroup SIR models with nonlinear incidence. Computers & Mathematics with Applications 2010, 60(8), 2286–2291. [Google Scholar]
  31. Koufi, A.L.; Adnani, J.; Bennar, A.; Yousfi, N. Analysis of a stochastic SIR model with vaccination and nonlinear incidence rate. Computers & Mathematics with Applications 2019, 9275051. [Google Scholar]
  32. Zhou, J.; Yang, Y.; Zhang, T. Global stability of a discrete multigroup SIR model with nonlinear incidence rate. Mathematical Methods in the Applied Sciences 2017, 40(14), 5370–5379. [Google Scholar] [CrossRef]
  33. Jin, Y.; Wang, W.; Xiao, S. An SIRS model with a nonlinear incidence rate. Chaos, Solitons and Fractals 2007, 34(5), 1482–1497. [Google Scholar] [CrossRef]
  34. Li, T.; Zhang, F.; Liu, H.; Chen, Y. Threshold dynamics of an SIRS model with nonlinear incidence rate and transfer from infectious to susceptible. Applied Mathematics Letters 2017, 70, 52–57. [Google Scholar] [CrossRef]
  35. Nudee, K.; Chinviriyasit, S.; Chinviriyasit, W. The effect of backward bifurcation in controlling measles transmission by vaccination. Chaos, Solitons and Fractals 2019, 123, 400–412. [Google Scholar] [CrossRef]
  36. Shan, C.; Shu, H. Bifurcations and complex dynamics of an SIR model with the impact of the number of hospital beds. Journal of Differential Equations 2014, 257, 1662–1688. [Google Scholar] [CrossRef]
  37. Turkyilmazoglu, M. A restricted epidemic SIR model with elementary solutions. Physica A 2022, 600, 127570. [Google Scholar] [CrossRef]
  38. Liu, L.; Jiang, D.; Hayat, T. Dynamics of an SIR epidemic model with varying population sizes and regime switching in a two patch setting. Physica A 2021, 574, 125992. [Google Scholar] [CrossRef]
  39. Arazi, R.; Feigel, A. Discontinuous transitions of social distancing in the SIR model. Physica A 2021, 566, 125632. [Google Scholar] [CrossRef]
  40. d’Onofrio, A.; Manfredi, P. Behavioral SIR models with incidence-based social-distancing. Chaos, Solitons and Fractals 2022, 159, 112072. [Google Scholar] [CrossRef]
  41. Giménez-Mujica, U.J.; Anzo-Hernández, A.; Velázquez-Castro, J. Epidemic local final size in a metapopulation network as indicator of geographical priority for control strategies in SIR type diseases. Mathematical Biosciences 2022, 343, 108730. [Google Scholar] [CrossRef]
  42. Yang, H.; Wang, Y.; Kundu, S.; Song, Z.; Zhang, Z. Dynamics of an SIR epidemic model incorporating time delay and convex incidence rate. Results in Physics 2022, 32, 105025. [Google Scholar] [CrossRef]
  43. Ndii, M. Z.; Beay, L. K.; Anggriani, N.; Nukul, K. N.; Zhang, Djahi, B.S. Estimating the Time Reproduction Number in Kupang City Indonesia, 2016–2020, and Assessing the Effects of Vaccination and Different Wolbachia Strains on Dengue Transmission Dynamics. Mathematics 2022, 10(12), 2075. [Google Scholar] [CrossRef]
Preprints 83664 g001
Figure 1. Scheme of disease transmission.
Figure 1. Scheme of disease transmission.
Preprints 83664 g001
Preprints 83664 g002
Figure 2. Stability of co-existing point E * with initial values (150,60,0,0).
Figure 2. Stability of co-existing point E * with initial values (150,60,0,0).
Preprints 83664 g002
Preprints 83664 g003
Figure 3. Stability of non-endemic point E 0 with initial values (150,60,0,0).
Figure 3. Stability of non-endemic point E 0 with initial values (150,60,0,0).
Preprints 83664 g003
Preprints 83664 g004
Figure 4. Bifurcation diagram for the force of infection β 1 of the model (1).
Figure 4. Bifurcation diagram for the force of infection β 1 of the model (1).
Preprints 83664 g004
Preprints 83664 g005
Figure 5. Simulation the effects of k, β 1 , α 0 , and α 1 on R and W
Figure 5. Simulation the effects of k, β 1 , α 0 , and α 1 on R and W
Preprints 83664 g005
Preprints 83664 g006
Figure 6. Effects of ξ on R and W
Figure 6. Effects of ξ on R and W
Preprints 83664 g006
Preprints 83664 g007
Figure 7. Control function.
Figure 7. Control function.
Preprints 83664 g007
Preprints 83664 g008
Figure 8. Simulation the impact of optimal control. Note: Red represents the simulation without control, and blue represents the simulation with control.
Figure 8. Simulation the impact of optimal control. Note: Red represents the simulation without control, and blue represents the simulation with control.
Preprints 83664 g008
Table 2. Parameters Values
Table 2. Parameters Values
Parameter Value
b 2.3
β 1 0.01
β 2 0.01
μ 0.005
γ 0.2
ξ 0.02
α 0 0.2
α 1 0.21
k 2
m 0.2
Table 3. Initial Values of Each Compartment
Table 3. Initial Values of Each Compartment
Compartment S ( 0 ) I ( 0 ) R ( 0 ) W ( 0 )
Initial Values 150 60 0 0
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
Copyright: This open access article is published under a Creative Commons CC BY 4.0 license, which permit the free download, distribution, and reuse, provided that the author and preprint are cited in any reuse.
Prerpints.org logo

Preprints.org is a free preprint server supported by MDPI in Basel, Switzerland.

facebook logotwitter logolinkedin logo
MDPI Initiatives
Important Links
Subscribe

Disclaimer

Terms of Use

Privacy Policy

© 2024 MDPI (Basel, Switzerland) unless otherwise stated