
A generalized "SVEIR" epidemic model with general nonlinear incidence rate has been proposed as a candidate model for measles virus dynamics. The basic reproduction number R, an important epidemiologic index, was calculated using the next generation matrix method. The existence and uniqueness of the steady states, namely, disease-free equilibrium (E0) and endemic equilibrium (E1) was studied. Therefore, the local and global stability analysis are carried out. It is proved that E0 is locally asymptotically stable once R is less than. However, if R>1 then E0 is unstable. We proved also that E1 is locally asymptotically stable once R>1. The global stability of both equilibrium E0 and E1 is discussed where we proved that E0 is globally asymptotically stable once R≤1, and E1 is globally asymptotically stable once R>1. The sensitivity analysis of the basic reproduction number R with respect to the model parameters is carried out. In a second step, a vaccination strategy related to this model will be considered to optimise the infected and exposed individuals. We formulated a nonlinear optimal control problem and the existence, uniqueness and the characterisation of the optimal solution was discussed. An algorithm inspired from the Gauss-Seidel method was used to resolve the optimal control problem. Some numerical tests was given confirming the obtained theoretical results.
Citation: Miled El Hajji, Amer Hassan Albargi. A mathematical investigation of an 'SVEIR' epidemic model for the measles transmission[J]. Mathematical Biosciences and Engineering, 2022, 19(3): 2853-2875. doi: 10.3934/mbe.2022131
[1] | Ke Guo, Wanbiao Ma . Global dynamics of an SI epidemic model with nonlinear incidence rate, feedback controls and time delays. Mathematical Biosciences and Engineering, 2021, 18(1): 643-672. doi: 10.3934/mbe.2021035 |
[2] | Yu Ji . Global stability of a multiple delayed viral infection model with general incidence rate and an application to HIV infection. Mathematical Biosciences and Engineering, 2015, 12(3): 525-536. doi: 10.3934/mbe.2015.12.525 |
[3] | Jinliang Wang, Gang Huang, Yasuhiro Takeuchi, Shengqiang Liu . Sveir epidemiological model with varying infectivity and distributed delays. Mathematical Biosciences and Engineering, 2011, 8(3): 875-888. doi: 10.3934/mbe.2011.8.875 |
[4] | Pengyan Liu, Hong-Xu Li . Global behavior of a multi-group SEIR epidemic model with age structure and spatial diffusion. Mathematical Biosciences and Engineering, 2020, 17(6): 7248-7273. doi: 10.3934/mbe.2020372 |
[5] | Pengfei Liu, Yantao Luo, Zhidong Teng . Role of media coverage in a SVEIR-I epidemic model with nonlinear incidence and spatial heterogeneous environment. Mathematical Biosciences and Engineering, 2023, 20(9): 15641-15671. doi: 10.3934/mbe.2023698 |
[6] | Maoxing Liu, Yuhang Li . Dynamics analysis of an SVEIR epidemic model in a patchy environment. Mathematical Biosciences and Engineering, 2023, 20(9): 16962-16977. doi: 10.3934/mbe.2023756 |
[7] | Tingting Xue, Long Zhang, Xiaolin Fan . Dynamic modeling and analysis of Hepatitis B epidemic with general incidence. Mathematical Biosciences and Engineering, 2023, 20(6): 10883-10908. doi: 10.3934/mbe.2023483 |
[8] | N. H. AlShamrani, A. M. Elaiw . Stability of an adaptive immunity viral infection model with multi-stages of infected cells and two routes of infection. Mathematical Biosciences and Engineering, 2020, 17(1): 575-605. doi: 10.3934/mbe.2020030 |
[9] | Thomas Torku, Abdul Khaliq, Fathalla Rihan . SEINN: A deep learning algorithm for the stochastic epidemic model. Mathematical Biosciences and Engineering, 2023, 20(9): 16330-16361. doi: 10.3934/mbe.2023729 |
[10] | Xiao-Min Huang, Xiang-ShengWang . Traveling waves of di usive disease models with time delay and degeneracy. Mathematical Biosciences and Engineering, 2019, 16(4): 2391-2410. doi: 10.3934/mbe.2019120 |
A generalized "SVEIR" epidemic model with general nonlinear incidence rate has been proposed as a candidate model for measles virus dynamics. The basic reproduction number R, an important epidemiologic index, was calculated using the next generation matrix method. The existence and uniqueness of the steady states, namely, disease-free equilibrium (E0) and endemic equilibrium (E1) was studied. Therefore, the local and global stability analysis are carried out. It is proved that E0 is locally asymptotically stable once R is less than. However, if R>1 then E0 is unstable. We proved also that E1 is locally asymptotically stable once R>1. The global stability of both equilibrium E0 and E1 is discussed where we proved that E0 is globally asymptotically stable once R≤1, and E1 is globally asymptotically stable once R>1. The sensitivity analysis of the basic reproduction number R with respect to the model parameters is carried out. In a second step, a vaccination strategy related to this model will be considered to optimise the infected and exposed individuals. We formulated a nonlinear optimal control problem and the existence, uniqueness and the characterisation of the optimal solution was discussed. An algorithm inspired from the Gauss-Seidel method was used to resolve the optimal control problem. Some numerical tests was given confirming the obtained theoretical results.
Measles is an infectious disease caused by a virus of the genus Morbilivirus and belonging to a group RNA viruses. The extremely contagious measles virus is spread when people cough or sneeze, through close contact between people or through direct contact with nasal, laryngeal secretions or with contaminated objects [1,2]. Symptoms of measles are characterized by a high fever accompanied by cough, runny nose, red and watery eyes, and Koplik's spot, which appear about 10–12 days after exposed to the virus [1]. In 2018, more than 140,000 people worldwide died from measles in part large are children under 5 years of age [3]. Measles immunization is given to infants aged 9 months because babies under that age still have innate immunity from their mother. Measles virus quickly killed by ultraviolet rays, chemicals, acids, and heating. Hence, early handling against measles can be done after going through laboratory confirmation by doing serological examination (patient blood collection/blood serum) or virological examination (patient urine collection) [2].
In the current context where we live serious health crises both in the industrialized countries and those who are less industrialized, one of the disciplinary fields of predilection embraced by mathematical modeling is public health. Several research teams are interested in the potential of mathematical models in the development of new tools and methods to control the spread of a disease [4]. Mathematical modeling of infectious diseases makes it possible to study the unexplained questions of an epidemic in a population, its dynamics of transmission, the probability that the epidemic will spread or die out in the susceptible population, the best vaccination strategy that would be effective to be undertaken by governments and any possibility of treatment of the disease. Develop, validate and use mathematical models to analyze and predict the dynamics of infectious diseases such as HIV, tuberculosis and hepatitis using the available data [5]. In the context of predictions and insights for the time-evolution related to infectious diseases, several researchers developed mathematical compartmental models; for example, see [6,7,8,9,10,11,12,13,14,15].
The spread of infectious diseases can be modeled with an epidemiological model, using either deterministic or stochastic approach [16,17,18,19,20,21]. In the health sector, epidemiological models can be used for proposing and testing theories, as well as for comparing, planning, implementing, and evaluating various detection, prevention, or control programs [22]. One of the basic models for the spread of disease was proposed by Kermack and McKendrick [23] in 1927. In their original model, Kermack and McKendrick introduced a model of dividing differential equations population into three compartments, namely the infected population compartment (I), the susceptible population compartment (S), and the recovered population compartment (R), known as the Susceptible (susceptible)-Infected (infected)-Recovered ("SIR") model [23]. The "SIR" model is used to model disease assuming that it is an individual have permanent immunity or temporary immunity for a very long period of time. Measles epidemic models have been studied by several researchers using deterministic models, as in [2,7,24,25]. In [25], Momoh et al. used the "SEIR" (Susceptible, Exposed, Infected, Recovered) to analyze the effect of testing and therapy for measles during the incubation period (exposed) against the dynamics of the spread of measles. Furthermore, Edward et al. [24] using the "SVEIR" model (Susceptible, Vaccinated, Exposed, Infected, Recovered) taking into account the proportion of immigrants who have been vaccinated, the proportion of the population that has been vaccinated twice, and the rate of effectiveness of the first vaccine. It is modified by Aldila and Asrianti [7] by differentiating the compartments of the population that have been vaccinated once and the population compartment that has been vaccinated twice, as well as the population compartment that has been vaccinated do quarantine. Fakhruddin et al. [2] modified the "SIR" model by adding a hospitalized (Hospitalized) infected population compartments.
The application of the "SVEIR" model to the spread of measles as reviewed in this article is in line with government efforts to control, reduce, and even eliminate the spread of measles human population, by carrying out a mass measles vaccine since 2000 [2]. Based on Therefore, the model reviewed in this article modifies the model proposed by Momoh et al. [25] with add vaccinated population compartments, and analyze the dynamics of the spread measles by looking at the influence of the proportion of susceptible population who have been vaccinated and the rate of infection when susceptible individuals interact with infected individuals. In [26], Wei et al. proposed an "SVEIR" epidemic model where a pulse vaccination was given for each individual. The strategy is to vaccinate each individual several times by separating doses of a definite time of a certain age group. Therefore, in the present article, we shall revisit and analyse the model first proposed by Hethcote [18], then studied by Gumel et al. [12], Wei et al. [26], Adda et al. [27], and then by Nkamba et al. [28] but by considering general nonlinear increasing incidence rate with respect to the infected individuals. Our choice is motivated by the fact that the number of effective contacts between infected individuals and susceptible individuals or between vaccinated individuals and infected individuals may increase at high infective levels due to crowding of infected individuals. Further, we proposed an optimal vaccination strategy relative to this epidemic model in order to minimize the number of both infected and exposed individuals.
The organization of the paper is as follows: In a first step, we present, in Section 2, a description of the investigated model. Further, the basic properties of the model for the existence of the solutions are obtained. The basic reproduction number R was calculated using the next generation matrix method. Then existence, uniqueness of the model equilibria was studied. Local stability of the equilibria of the proposed model was investigated in Section 3. The global stability and asymptotic behaviour of the solution are discussed in Section 4. In Section 5, the sensitivity analysis of the basic reproduction number R with respect to the model parameters is carried out. In Section 6, we proposed and formulated an optimal vaccination strategy in order to optimize the number of infected and exposed individuals. We discussed the existence, uniqueness and characterisation of the optimal solution. Finally, in Section 7, we give some numerical simulations, validating the theoretical results where we applied an appropriated numerical scheme inspired from the competitive Gauss-Seidel like implicit difference method.
We assume that the total population is divided into five classes (compartments) according to their status concerning disease: susceptible S(t), vaccinated V(t), exposed E(t), infected I(t), and recovered R(t) classes. The susceptible class, S(t), represents the individuals that having the risk to catch the infection because of the close contact with the infected individuals. The vaccinated class, V(t), consists of those individuals who may be given an imperfect vaccine that reduces their susceptibility to the disease; still, they can get the infection under some suitable conditions. The infected class, I(t), represents the individuals who already cached the disease, and they are able to infect the susceptible individuals for some suitable conditions. The recovered class, R(t), represents the individuals who are already infected, and now they recovered their healthy due to either treatment or immunization. Schematically, the spread of measles using the "SVEIR" model is presented in the following transfer diagram.
To build a deterministic model for the "SVEIR" spread of measles, several things are assumed following.
● Λ is the recruitment rate of susceptible individuals.
● p describe the rate of vaccination.
● θ∈[0,1] is the effectiveness rate of the vaccination.
● Each compartment has a specific natural mortality rate. ms, mv, me, mi and mr are the mortality rates of susceptible, vaccinated, exposed, infected and recovered individuals, respectively.
● Measles can cause death (fatal) with a mortality rate of ζ.
● f represents the saturated incidence rate.
● ε the rate at which an exposed individuals become infectious. 1/ε is the average latency time spent in compartment E.
● γ the rate at which infectious agents recover their health. 1/γ is the average duration elapsed in compartment I.
● Individuals who have been vaccinated will no gain permanent immunity so they will can be infected again.
● Measles has an incubation period so that when the incubation period ends, individuals in the Exposed (E) compartment will move to the Infected (I) compartment at a rate of ε.
● The recovery rate of individuals infected with measles is γ.
● Recovered individuals will not re-infect.
The generalized "SVEIR" mathematical model of the ones given in [12,26,27,28] for the spread of the epidemic based on the transfer diagram in Figure 1 is given by:
{˙S=Λ−(ms+p)S−f(I)S,˙V=pS−mvV−θf(I)V,˙E=f(I)(S+θV)−(me+ε)E,˙I=εE−(mi+ζ+γ)I,˙R=γI−mrR, | (2.1) |
with positive initial condition (S(0),V(0),E(0),I(0),R(0))∈R5+.
Variable | Description | Parameter | Description |
S(t) | Susceptible individuals | f | Saturated incidence rate |
V(t) | Vaccinated individuals | Λ | Susceptible recruitment rate |
E(t) | Exposed individuals | ε | Infection transmission rate |
I(t) | Infected individuals | γ | Recovery rate |
R(t) | Recovered individuals | θ∈[0,1] | Effectiveness rate of the vaccination |
Parameter | Description | Parameter | Description |
mv | Mortality rate of vaccinated individuals | ms | mortality rate of susceptible individuals |
me | Mortality rate of exposed individuals | p | vaccination rate |
mi+ζ | Total mortality rate of infected individuals | 1ε | Average latency time spent in compartment E |
mr | Mortality rate of recovered individuals | 1γ | Average duration elapsed in compartment I |
Assumption 1. The function f is non-negative C1(R+), increasing bounded and concave such that f(0)=0.
Any function f satisfying Assumption 1 satisfies the following lemma.
Lemma 1. f′(I)I<f(I)≤f′(0)I for all I>0.
Proof. For I∈R+, let h1(I)=f(I)−If′(I). Since f is an increasing and concave function then f′(I)≥0 and f″(I)≤0. Therefore h′1(I)=−If″(I)≥0 and h1(I)≥h1(0)=0 or also f(I)≥If′(I). Similarly, let h2(I)=f(I)−If′(0) then h′2(I)=f′(I)−f′(0)≤0 since f is concave. Then h2(I)≤h2(0)=0 then f(I)≤If′(0).
Remark 1. The Monod (or Holling's type Ⅱ) functions satisfy the assumption 1 and then can express the transmission rate.
f(I)=ˉfIk+I |
ˉf designs the maximum transmission rate and k represents the Monod constant.
In order to prove that the system (2.1) is well-posed, we will give some general properties. Let m=min(ms,mv,me,mi+ζ,mr), then
Proposition 1. The compact set Ω1={(S,V,E,I,R)∈R5+/S+V+E+I+R≤Λm} is positively invariant for system (2.1).
Proof. Since ˙S|S=0=Λ>0, ˙V|V=0=pS>0, ˙E|E=0=f(I)(S+θV)>0, ˙I|I=0=εE>0 and ˙R|R=0=γI>0, then the solution of system (2.1) is non-negative.
Consider a new variable T=S+V+E+I+R−Λm, then by adding all equations of (2.1), T satisfies:
˙T=˙S+˙V+˙E+˙I+˙R=Λ−msS−mvV−meE−(mi+ζ)I−mrR≤Λ−m(S+V+E+I+R)≤m(Λm−S−V−E−I−R)≤−mT. |
Hence
T(t)≤T(0)e−mt. | (2.2) |
Since all variables are non-negative then all variables are bounded and thus Ω1 is positively invariant for system (2.1).
In epidemiology, a major health index calculating the average number of people an infectious person infects while they are contagious, is called the basic reproduction number, and it is denoted R. Recall that R is the number of secondary infections produced by a single typical infection in a susceptible population, a very common question is how to define R when there are several types of infected individuals, for example we can have an atypical infection in a disease vector-borne such as malaria, or sexually transmitted infections when there are large asymmetries in transmission, such as HIV, or a multiple-host pathogen such as influenza. The key concept is that we now need to average the expected number of news infections in our case. If there are several compartments representing infectious individuals, the next-generation matrix method introduced by Diekmann [4] and developed later by van den Driessche and Watmough [5] is used to calculate the basic reproduction number R. Let S0=Λms+p,V0=pΛmv(ms+p),E0=0,ˉI0=0 and R0=0. The Jacobian of the sub-system (E,I) of the original system (2.1) at the disease-free equilibrium point E0=(S0,V0,E0,I0,R0) is given by
J(E0)=((me+ε)f′(0)(S0+θV0)−ε(mi+ζ+γ))=F+V |
where F=(0f′(0)(S0+θV0)00) and V=((me+ε)0−ε(mi+ζ+γ)).
The determinant of V is given by det(V)=(me+ε)(mi+ζ+γ)>0 and thus the inverse matrix of V is
V−1=(1(mi+ζ+γ)0ε(me+ε)(mi+ζ+γ)1(me+ε)) |
and the next-generation matrix is given by
FV−1=(εf′(0)(S0+θV0)(me+ε)(mi+ζ+γ)f′(0)(S0+θV0)(me+ε)00). |
Then, the basic reproduction number of system (2.1) is calculated as the spectral radius of the matrix FV−1:
R=εf′(0)(S0+θV0)(me+ε)(mi+ζ+γ)=Λεf′(0)(mv+θp)mv(me+ε)(mi+ζ+γ)(ms+p). | (2.3) |
Note that in absence of vaccination (p=0), the system (2.1) becomes the standard SEIR model with
R0=R|p=0=Λεf′(0)ms(me+ε)(mi+ζ+γ). | (2.4) |
Proposition 2. System (2.1) admits a trivial (disease-free) equilibrium point E0=(S0,V0,E0,I0,R0) and if R>1, system (2.1) admits a unique non-trivial (endemic) equilibrium point E1=(S1,V1,E1,I1,R1) with S1,V1,E1,I1,R1>0.
Proof. An equilibrium point for the system (2.1) satisfies
{0=Λ−(ms+p)S−f(I)S,0=pS−mvV−θf(I)V,0=f(I)(S+θV)−(me+ε)E,0=εE−(mi+ζ+γ)I,0=γI−mrR, | (2.5) |
which reduces to
{S=Λ(ms+p+f(I)),V=pS(mv+θf(I))=pΛ(mv+θf(I))(ms+p+f(I)),E=f(I)(me+ε)(S+θV)=Λf(I)(mv+θp+θf(I))(me+ε)(ms+p+f(I))(mv+θf(I)),I=ε(mi+ζ+γ)E=εΛf(I)(mv+θp+θf(I))(mi+ζ+γ)(me+ε)(ms+p+f(I))(mv+θf(I)),R=γmrI=γεΛf(I)(mv+θp+θf(I))mr(mi+ζ+γ)(me+ε)(ms+p+f(I))(mv+θf(I)). | (2.6) |
From the fourth equation of (2.6) one deduces that
I=εΛf(I)(mv+θp+θf(I))(mi+ζ+γ)(me+ε)(ms+p+f(I))(mv+θf(I)). |
Since all parameters are non-negative then either I=0 or
εΛf(I)(mv+θp+θf(I))=I((mi+γ)(me+ε)(ms+p+f(I))(mv+θf(I))). |
● If I=0 then R=0,E=0,S=Λms+p=ˉS and V=pSmv=pΛmv(ms+p)=ˉV. This equilibrium known as the disease-free equilibrium denoted here by E0=(S0,V0,E0,I0,R0).
● For I≠0, define the continuous function g
g(I)=1−εΛf(I)I(mv+θp+θf(I))(mi+ζ+γ)(me+ε)(ms+p+f(I))(mv+θf(I)). |
The derivative of g is given by
g′(I)=−εΛ(mi+ζ+γ)(me+ε)(ms+p+f(I))2(mv+θf(I))2×([f′(I)I−f(I)I2(mv+θp+θf(I))+θf′(I)f(I)I](ms+p+f(I))(mv+θf(I))−f′(I)f(I)I(mv+θp+θf(I))(mv+θf(I)+θ(ms+p+f(I))))=−εΛf′(I)I−f(I)I2(mv+θp+θf(I))(mi+ζ+γ)(me+ε)(ms+p+f(I))(mv+θf(I))−εΛf′(I)f(I)I(mi+ζ+γ)(me+ε)(ms+p+f(I))2(mv+θf(I))2×(θ(ms+p+f(I))(mv+θf(I))−(mv+θp+θf(I))(mv+θf(I)+θ(ms+p+f(I))))=−εΛf′(I)I−f(I)I2(mv+θp+θf(I))(mi+ζ+γ)(me+ε)(ms+p+f(I))(mv+θf(I))+εΛf′(I)f(I)I((mv+θp+θf(I))(mv+θf(I))+θ2p(ms+p+f(I))))(mi+ζ+γ)(me+ε)(ms+p+f(I))2(mv+θf(I))2. |
Since all parameters are non-negative and since the function f satisfies f′(I)I<f(I) for all I>0, one can easily deduce that the function g is an increasing function.
A simple calculus gives
limI↦0g(I)=1−εΛlimI↦0f(I)I(mv+θp+θlimI↦0f(I))(mi+ζ+γ)(me+ε)(ms+p+limI↦0f(I))(mv+θlimI↦0f(I))=1−εΛlimI↦0f(I)I(mv+θp)mv(mi+ζ+γ)(me+ε)(ms+p)=1−εΛf′(0)(mv+θpmv(mi+ζ+γ)(me+ε)(ms+p)=1−R, |
and
g(Λm)=1−εmf(Λm)(mv+θp+θf(Λm))(mi+ζ+γ)(me+ε)(ms+p+f(Λm))(mv+θf(Λm))>1−εmf(Λm)(mv+θp+θf(Λm))(mi+ζ)ε(ms+p+f(Λm))(mv+θf(Λm))>1−εmf(Λm)(mv+θp+θf(Λm))mε(ms+p+f(Λm))(mv+θf(Λm))>1−εmf(Λm)(mv+θp+θf(Λm))mεf(Λm)(mv+θf(Λm))+mεpθf(Λm)+mεmsmv=εmsmvεf(Λm)(mv+θf(Λm))+εpθf(Λm)+εmsmv>0 |
Since R>1, limI↦0g(I)<0 and g(Λm)>0, the equation g(I)=0 admits a unique solution I1 in (0,Λm) and then the uniqueness of the endemic equilibrium E1=(S1,V1,E1,I1,R1).
In the section below, we discussed the local stability of equilibria with respect to the value of the basic reproduction number R.
Theorem 1. If R<1, then the disease-free equilibrium (DFE) point E0 is locally asymptotically stable and it is unstable if R>1.
Proof. The Jacobian matrix associated to the system (2.1) at a given point (S,V,E,I,R) is calculated as follows:
J=(−(ms+p)−f(I)00−f′(I)S0p−mv−θf(I)0−θf′(I)V0f(I)θf(I)−(me+ε)f′(I)(S+θV)000ε−(mi+ζ+γ)0000γ−mr) |
It's value, evaluated at E0, is given by:
ˉJ=(−ms−p00−f′(0)S00p−mv0−θf′(0)V0000−me−εf′(0)(S0+θV0)000ε−mi−ζ−γ0000γ−mr) |
which admits five distinguish eigenvalues given by λ1=−(ms+p)<0,λ2=−mv<0,λ3=−mr<0 and two other eigenvalues corresponding to the sub-matrix
(−(me+ε)f′(0)(S0+θV0)ε−(mi+ζ+γ)). |
The corresponding characteristic polynomial is calculated as follows
P(λ)=λ2+(me+mi+ζ+ε+γ)λ+(1−R)(me+ε)(mi+ζ+γ). |
Only if R<1, then the roots of P(.) have negative real parts. Then E0 is locally asymptotically stable only if R<1. E0 is unstable once R>1.
In the next section, the global stability behavior of the equilibria will be discussed.
Corollary 1. The set Ω2={(S,V,E,I,R)∈R5+/S+V+E+I+R≤Λm;S≤S0,V≤V0} is positively invariant for system (2.1).
Proof. It is already proved that Ω1 is invariant for system (2.1). Note that ˙S(t)<0 for S(t)>S0 therefore lim infS(t)≤S0. Let ξ>0 be a constant and S(0) be an initial condition. Then ∃T≥0;S(t)≤S0+ξ∀t≥T. Therefore
˙V(t)<p(S0+ξ)−mvVfor allt≥T. |
Consider ˉV=V−p(S0+ξ)mv which satisfies
˙ˉV(t)<−mvˉVfor allt≥T. | (4.1) |
Then
˙ˉV(t)<˜V(0)e−mvtfor allt≥T. | (4.2) |
Therefore lim infˉV(t)≤0 and lim infV(t)≤p(S0+ξ)mv for all ξ>0 thus lim infV(t)≤pS0mv=V0. This completes the proof.
Theorem 2. E0 is globally asymptotically stable if R≤1 however it is unstable if R>1.
Proof. Let the function F0 given by:
F0=εE+(me+ε)I. |
Since f(I)≤f′(0)I, the time-derivative of F0 is given by:
˙F0=ε˙E+(me+ε)˙I=ε(f(I)(S+θV)−(me+ε)E)+(me+ε)(εE−(mi+ζ+γ)I)=εf(I)(S+θV)−(me+ε)(mi+γ)I≤εf′(0)I(S+θV)−(me+ε)(mi+ζ+γ)I≤(εf′(0)(S+θV)−(me+ε)(mi+ζ+γ))I≤(me+ε)(mi+ζ+γ)(εf′(0)(S+θV)(me+ε)(mi+ζ+γ)−1)I≤(me+ε)(mi+ζ+γ)(εf′(0)(ˉS+θˉV)(me+ε)(mi+ζ+γ)−1)I,sinceS≤ˉS,V≤ˉV=(me+ε)(mi+ζ+γ)(R−1)I,∀(S,V,E,I,R)∈Ω2. |
If R≤1, then ˙F0≤0 for all S,V,E,I,R>0. Let W0={(S,V,E,I,R):˙F0=0}={E0}. Then by LaSalle's invariance principle [29] E0 is globally asymptotically stable once R≤1 (see [9,20,21,30] for some applications). Then the solution of system (2.1) converges to E0 as t→+∞.
Hereafter, we study the global stability analysis for the endemic equilibrium E1.
Theorem 3. If R>1, then the endemic equilibrium E1 is globally asymptotically stable.
Proof. Let the function F1, given by
F1=(S−S1ln(SS1))+(V−V1ln(VV1))+(E−E1ln(EE1))+me+εε(I−I1ln(II1)). |
The function F1 admits its minimum value F1min=S1+V1+E1+me+εεI1 once S=S1,V=V1,E=E1,I=I1. The time-derivative of F1, along solutions of system (2.1) is given by
˙F1=(1−S1S)˙S+(1−V1V)˙V+(1−E1E)˙E+me+εε(1−I1I)˙I=(1−S1S)(Λ−(ms+p)S−f(I)S)+(1−V1V)(pS−mvV−θf(I)V)+(1−E1E)(f(I)(S+θV)−(me+ε)E)+me+εε(1−I1I)(εE−(mi+ζ+γ)I)=Λ−msS−ΛS1S+msS1+pS1+f(I)S1−mvV−pSV1V+mvV1+θf(I)V1−f(I)(S+θV)E1E+(me+ε)E1−me+εε(mi+ζ+γ)I−(me+ε)EI1I+me+εε(mi+ζ+γ)I1=msS1−msS−msS∗2S+msS1+mvV1−mvS1SV1+mvV1−mvV+mvV1+mvV1−mvV∗2VSS1+f(I1)(S1+θV1)−f(I1)S∗2S−θf(I1)V1S1S+θf(I1)V1+f(I1)S1−θf(I1)V∗2VSS1 |
+θf(I)V1−f(I)(S+θV)E1E+f(I)(S1+θV1)−me+εε(mi+ζ+γ)I−(me+ε)EI1I+me+εε(mi+ζ+γ)I1=msS1(2−SS1−S1S)+mvV1(3−S1S−VV1−V1VSS1)+3f(I1)S1−f(I1)S∗2S−f(I1)S1EE1I1I−f(I)SE1E−3θf(I1)V1+θf(I1)V1−θf(I1)V∗2VSS1−θf(I)VE1E−θf(I1)V1EE1I1I−θf(I1)V1S1S−me+εε(mi+ζ+γ)I+f(I)(S1+θV1)=msS1(2−SS1−S1S)+mvV1(3−S1S−VV1−V1VSS1)+f(I1)S1(3−S1S−II1SS1E1E−EE1I1I)+θf(I1)V1(4−V1VSS1−II1VV1E1E−EE1I1I−S1S). |
Since
x1+x2+x3+⋯+xn≥nn√x1.x2.x3⋯xn,x1,x2,x3,⋯,xn≥0 | (4.3) |
then (2−SS1−S1S)≤0,(3−S1S−VV1−V1VSS1)≤0, (3−S1S−II1SS1E1E−EE1I1I)≤0, (4−V1VSS1−II1VV1E1E−EE1I1I−S1S)≤0. Therefore ˙F2≤0. Thank's to Lyapunov theorem, E1 is stable.
Now, we have to show the asymptotic stability of E1 using the Lasalle invariance principle.
Define A=(2−SS1−S1S),B=(3−S1S−VV1−V1VSS1), C=(3−S1S−II1SS1E1E−EE1I1I),
D=(4−V1VSS1−II1VV1E1E−EE1I1I−S1S). Then ˙F1(S,V,E,I,R)=0 if and only if A=B=C=D=0.
A=0 means that S=S1 and B=0 means that V=V1. If moreover C=0 then EE1=II1. To conclude
˙F1(S,V,E,I,R)=0⇔S=S1,V=V1,EE1=II1. |
Let a=EE1=II1, then E=aE1 and I=aI1.
The endemic equilibrium satisfies
|Λ=(ms+p)S1+f(I1)S1,pS1=mvV1+θf(I1)V1,f(I1)(S1+θV1)=(me+ε)E1,εE1=(mi+ζ+γ)I1. |
Then
pS1=mvV1+θf(I1)V1,a=1 |
and thus I=I1 and E=E1. Finally
˙F1(S,V,E,I,R)=0⇔(S=S1,V=V1,E=E1,I=I1,R=R1). |
Therefore {(S,V,E,I,R)|˙F1=0}={E1}. Then by LaSalle's invariance principle [29] E1 is globally asymptotically stable once R>1 (see [9,20,21,30] for some applications).
The sensitivity analysis helps to determine how the variability of the model outputs can be attributed to different sources of variation of its input parameters [31,32]. It can thus have various objectives:
● Identify the input parameters that most contribute to an output of interest from the model. Consequently, this makes it possible to prioritize the parameters to be estimated rather than to estimate all the uncertain parameters in order to save simulation time considerable (most common goal).
● Detect and quantify the effects of interactions between input parameters.
● Determine possible simplifications of the model.
In our case, the sensitivity analysis can tells us how important each parameter with respect to the disease transmission. Thus we have to discover the parameters that have the most high impact on the basic reproduction number R. These parameters will be targeted by any possible intervention strategies. A key parameter known as the normalized forward sensitivity index describe the ratio of the relative change in the variable to the relative change in this parameter. If the result is positive, this means that an increase in the parameter value induces an increase in the basic reproduction number value. However, a negative sensitivity index means that the parameter and the basic reproduction number are inversely proportional [33].
Definition 1. (see [31,32]). The normalized forward sensitivity index of R that depends differentiably on a parameter ρ is defined by
YRρ=∂R∂ρ×ρ|R|. |
For instance, YRρ=1 implies an increase (decrease) of ρ by y% increases (decreases) R by the same percentage. On the other hand, YRρ=−1 indicates an increase (decrease) of ρ by y% decreases (increases) R by y%. In this case, ρ is called a highly sensitive parameter.
Proposition 3. The explicit expression of R is given by
R=Λεf′(0)(mv+θp)mv(me+ε)(mi+ζ+γ)(ms+p). |
The basic reproduction number R depends on some parameters, we calculate analytical expressions for its sensitivity to each of these parameters by calculating the normalized forward sensitivity index as the following:
Proof. This follows immediately from Definition 1.
From the Table 2, it is to be noted that the parameters Λ,ε and θ are positive and hence an increase of the values of parameters Λ,ε and θ leads to an increase of the value of R. Note that the dependence of R on p depend on the sign of (θms−mv), and the dependence of R on the remaining parameters are negative and hence an increase of the values of these parameters leads to a decrease of the value of R.
Parameter | Sensitivity index of R | Sign |
Λ | YRΛ=∂R∂Λ×ΛR=1 | +ve |
ε | YRε=∂R∂ε×εR=meme+ε | +ve |
θ | YRθ=∂R∂θ×θR=θpmv+θp | +ve |
p | YRp=∂R∂p×pR=(θms−mv)p(ms+p)(mv+θp) | |
mv | YRmv=∂R∂mv×mvR=−θpmv+θp | -ve |
ms | YRms=∂R∂ms×msR=−msms+p | -ve |
me | YRme=∂R∂me×meR=−meme+ε | -ve |
mi | YRmi=∂R∂mi×miR=−mimi+ζ+γ | -ve |
ζ | YRζ=∂R∂ζ×ζR=−ζmi+ζ+γ | -ve |
γ | YRγ=∂R∂γ×γR=−γmi+ζ+γ | -ve |
The optimal control strategy for epidemic models remains a wide open research area [20,34,35]. In this section, we proposed an optimal vaccination strategy as the limit of preventive strategies for which vaccination effort is minimized. Since the fifth equation doesn't affect the other equations, let consider the reduced system given by
{˙S=Λ−(ms+p)S−f(I)S,˙V=pS−mvV−θf(I)V,˙E=f(I)(S+θV)−(me+ε)E,˙I=εE−(mi+ζ+γ)I, | (6.1) |
where p=p(t) is a function that will be considered as control.
Note that the compact set Ω1={(S,V,E,I)∈R4+/S+V+E+I≤Λm} is positively invariant for system (6.1).
Assume moreover that f is globally Lipschitz with an upper bound ˉf=supI>0f(I) and a Lipschitz constant Lf. Consider the admissible space
Vad={p|p(t)is measurable,pmin≤p(t)<pmax,0≤t≤T} |
where pmin and pmax are the bounds of the control.
The optimal control problem considered here is
minp∈VadJ(p)=minp∈Vad∫T0(a2p2(t)+αE(t)+βI(t))dtsubject to (6.1). | (6.2) |
Our objective for the work is to minimize the number of exposed E(t) and infective I(t) in the population and the cost of vaccination by using the optimal vaccination rate p(t). The quantities α and β denote the weight constants of the exposed and infective human population.
For φ=(S,V,E,I)t, the system (6.1) takes a new form
˙φ=Aφ+G(φ)=F(φ) | (6.3) |
where A=(−ms−p000p−mv0000−me−ε000ε−mi−ζ−γ) and G(φ)=(Λ−f(I)S−θf(I)V(S+θV)f(I)0).
Proposition 4. The continuous function F is uniformly Lipschitz.
Proof. The continuous function G is uniformly Lipschitz since
|G(φ1)−G(φ2)|=2|f(I1)S1−f(I2)S2|+2θ|f(I1)V1−f(I2)V2|≤2|S1(f(I1)−f(I2))+f(I2)(S1−S2)|+2θ|V1(f(I1)−f(I2))+f(I2)(V1−V2)|≤2|Λm(f(I1)−f(I2))+ˉf(S1−S2)|+2θ|Λm(f(I1)−f(I2))+ˉf(V1−V2)|≤2ΛmLf|I1−I2|+2ˉf|S1−S2|+2θΛmLf|I1−I2|+2ˉf|V1−V2|≤2max(ΛmLf,θΛmLf,ˉf)(|S1−S2|+|V1−V2|+|I1−I2|)≤M|φ1−φ2| |
where M=2max(ΛmLf,θΛmLf,ˉf). Since
|Aφ1−Aφ2|≤‖A‖|φ1−φ2| | (6.4) |
where ‖.‖ is the matrix norm, then
|F(φ1)−F(φ2)|≤K|φ1−φ2| | (6.5) |
where K=max(M,‖A‖) and the continuous function G is therefore uniformly Lipschitz.
Once the continuous function F is uniformly Lipschitz thus the existence and uniqueness of the solution of system (6.3).
By using Pontryagin's Maximum Principle [20,34,35,36], necessary conditions for our optimal control and corresponding states can be derived. Define the Hamiltonian H for the optimal control problems (6.1) and (6.2)
H(S,V,E,I,p,λ1,λ2,λ3,λ4)=a2p2+αE+βI+λ1˙S+λ2˙V+λ3˙E+λ4˙I=a2p2+αE+βI+λ1(Λ−(ms+p)S−f(I)S)+λ2(pS−mvV−θf(I)V)+λ3(f(I)(S+θV)−(me+ε)E)+λ4(εE−(mi+ζ+γ)I). |
For a given optimal control p∗, there exist adjoint functions λ1,λ2,λ3 and λ4 corresponding to the states S,V,E and I such that:
{˙λ1=−∂H∂S=p(λ1−λ2)+f(I)(λ1−λ3)+msλ1,˙λ2=−∂H∂V=θf(I)(λ2−λ3)+mvλ2,˙λ3=−∂H∂E=−α+ε(λ3−λ4)+meλ3,˙λ4=−∂H∂I=−β+f′(I)S(λ1−λ3)+θf′(I)V(λ2−λ3)+(mi+ζ+γ)λ4. | (6.6) |
Initial conditions of the adjoint variables are given at the time t = T by \lambda_1(T) = 0, \lambda_2(T) = 0, \lambda_3(T) = 0, \lambda_4(T) = 0 .
Since the Hamiltonian is minimized with respect to p^* thus the singular case could occur if the partial derivative function,
\begin{equation} \dfrac{\partial H}{\partial p} = a p - \lambda_1 S + \lambda_2 S, \end{equation} | (6.7) |
is zero on a non-trivial interval of time. Since the optimal control must be between an upper and lower bounds such that: \dfrac{\partial H}{\partial p} < 0 or > 0 .
In our case, the singular expression is given by
p_{\rm singular} (t) = \dfrac{S}{a}(\lambda_1-\lambda_2) \mbox{ if } p_{\min} \leq \dfrac{S}{a}(\lambda_1-\lambda_2) \leq p_{\max}. |
For the transmission rate of infection, we used Monod function (Known also as Holling's type Ⅱ): f(I) = \dfrac{\bar f I}{k+I} . \bar f designs the maximum transmission rate and k is a non-negative constant. Note that f is globally Lipschitz function on \mathbb{R}_+ . For all numerical simulations presented hereafter, we used the parameter values given in Table 3.
Parameter | \bar f | k | m_s | m_v | m_e | m_i | \zeta | m_r | p | \theta |
Value | 5 | 3 | 0.2 | 0.1 | 0.3 | 0.25 | 0.15 | 0.2 | 1 | 0.5 |
We give first some numerical simulations confirming the global stability of the equilibria of system (2.1).
For the both two cases where \mathcal{R} < 1 , the solution of system (2.1) converges asymptotically to \mathcal{E}_0 = (\bar S, \bar V, 0, 0, 0) (Figure 2). This confirms the global stability of \mathcal{E}_0 when \mathcal{R} \leq 1 . Similarly, the both last cases where \mathcal{R} > 1 , the solution of system (2.1) converges asymptotically to \mathcal{E}_1 (Figure 3). This validate the global stability of \mathcal{E}_1 when \mathcal{R} > 1 .
Concerning the numerical resolution of the control problem, we used an improved Gauss-Seidel implicit scheme for the state variables and a first-order backward-difference scheme for the adjoint-state variables [37]. For the numerical simulations of the optimal control problems (6.1) and (6.2), we used the same values as for the direct problem system (2.1) given in the Table 3 with \Lambda = 5, \epsilon = 0.5, \gamma = 0.5 with a variable p such that p(0) = 1 and with bounds p^{\min} = 0 and p^{\max} = 10 . We plot in Figure 4 the behaviours of the state variables S(t), V(t), E(t) and I(t) for different values of \alpha and \beta .
As seen in the Figure 4, vaccinated individuals increases about 37\% however the exposed and infected individuals decrease about 33 and 46\% , respectively, of their number at steady-state in Figure 3 (right). As it can be noted that no considerable influence of the values of \alpha and \beta on the behaviour of the variables S(t), V(t), E(t) and I(t) (Figure 4).
In this paper, we have proposed and studied a generalized "SVEIR" epidemic model with general nonlinear incidence rate. We discussed the stability of the equilibrium points (the trivial equilibrium point \mathcal{E}_0 and the nontrivial equilibrium point \mathcal{E}_1 ). \mathcal{E}_0 is locally asymptotically stable once the basic reproduction number \mathcal{R}\leq 1 and unstable once \mathcal{R} > 1 . We studied the uniqueness and existence of the nontrivial equilibrium point \mathcal{E}_1 therefore we studied its stability with respect to the basic reproduction number \mathcal{R} . We proved that \mathcal{E}_1 is locally asymptotically stable when \mathcal{R} is greater than 1. Later, we studied the global stability of both equilibrium points of system (2.1). We proved that \mathcal{E}_0 is globally asymptotically stable once \mathcal{R}\leq 1 , and that \mathcal{E}_1 is globally asymptotically stable once \mathcal{R} > 1 . The sensitivity analysis of the basic reproduction number \mathcal{R} with respect to the model parameters was carried out. In a second step, we proposed an optimal vaccination strategy to optimise the number of infected and exposed individuals. We formulated this strategy in terms of a nonlinear optimal control problem. We discussed the existence, uniqueness and the characterisation of the optimal solution. Finally, some numerical simulations that support the obtained analytical results are given.
The author would like to thank the editors and the anonymous reviewers whose invaluable comments and suggestions have greatly improved this manuscript.
The authors declare no conflict of interest.
Let a subdivision of the time interval [0, T] such that [0, T] = \bigcup\limits_{n = 0}^{N-1} [t_n, t_{n+1}] , t_n = n dt , dt = T/N . Let S^n, V^n, E^n, I^n , \lambda^n_1, \lambda^n_2, \lambda^n_3, \lambda^n_4 and p^n be an approximation of S(t) , V(t) , E(t) , I(t) , \lambda_1(t) , \lambda_2(t) , \lambda_3(t) , \lambda_4(t) and the control p(t) at the time t_n . S^0, V^0, E^0, I^0 , \lambda^0_1, \lambda^0_2, \lambda^0_3, \lambda^0_4 and p^0 as the state variables, the adjoint-state variables and the control at initial time. S^N, V^N, E^N, I^N , \lambda^N_1, \lambda^N_2, \lambda^N_3, \lambda^N_4 and p^N as the state variables, the adjoint-state variables and the control at final time T . The following appropriated scheme was then applied:
\begin{eqnarray*} \left\{ \begin{array}{ll} \dfrac{S^{n+1}-S^n}{\delta t} = \Lambda -(m_s+p^n)S^{n+1}-f(I^n)S^{n+1}, \\ \dfrac{V^{n+1}-V^n}{\delta t} = p^nS^{n+1}-m_vV^{n+1}-\theta f(I^n)V^{n+1}, \\ \dfrac{E^{n+1}-E^n}{\delta t} = f(I^n)(S^{n+1}+\theta V^{n+1})-(m_e+\varepsilon)E^{n+1}, \\ \dfrac{I^{n+1}-I^n}{\delta t} = \varepsilon E^{n+1} -(m_i+\zeta+\gamma)I^{n+1}, \\ \dfrac{\lambda_1^{N-n}-\lambda_1^{N-n-1}}{\delta t} = p^{n}(\lambda_1^{N-n-1}-\lambda_2^{N-n})+ f(I^{n+1})(\lambda_1^{N-n-1}-\lambda_3^{N-n}) +m_s \lambda_1^{N-n-1}, \\ \dfrac{\lambda_2^{N-n}-\lambda_2^{N-n-1}}{\delta t} = \theta f(I^{n+1})(\lambda_2^{N-n-1}-\lambda_3^{N-n})+m_v \lambda_2^{N-n-1}, \\ \dfrac{\lambda_3^{N-n}-\lambda_3^{N-n-1}}{\delta t} = -\alpha+\varepsilon(\lambda_3^{N-n-1}-\lambda_4^{N-n})+m_e \lambda_3^{N-n-1}, \\ \dfrac{\lambda_4^{N-n}-\lambda_4^{N-n-1}}{\delta t} = -\beta+f'(I^{n+1}) S^{n+1} (\lambda_1^{N-n-1}-\lambda_3^{N-n-1})\\ \qquad \qquad \qquad \qquad + \theta f'(I^{n+1})V^{n+1}(\lambda_2^{N-n-1}- \lambda_3^{N-n-1}) + (m_i+\zeta+\gamma)\lambda_4^{N-n-1} \end{array} \right. \label{OptimalSystem} \end{eqnarray*} |
Hence, the following algorithm was applied.
Algorithm 1: Optimal control resolution |
S^0 \gets S_0, V^0 \gets V_0, E^0 \gets E_0 , I^0 \gets I_0 , \lambda_1^N \gets 0, \lambda_2^N \gets 0, \lambda_3^N \gets 0, \lambda_4^N \gets 0 , for n = 0 to N-1 do \begin{eqnarray*} \left[ \begin{array}{llll} S^{n+1} \gets \dfrac{S^n+\delta t \Lambda}{1+\delta t (m_s+p^n)+ \delta tf(I^n)}, \\ V^{n+1} \gets \dfrac{V^n+\delta t p^nS^{n+1}}{1+\delta t m_v+ \delta t\theta f(I^n)}, \\ E^{n+1} \gets \dfrac{E^n+\delta t f(I^n)(S^{n+1}+\theta V^{n+1})}{1+ \delta t (m_e+\varepsilon)}, \\ I^{n+1} \gets \dfrac{I^n+\delta t \varepsilon E^{n+1}}{1+\delta t (m_i+\zeta+\gamma)}, \\ \lambda_1^{N-n-1} \gets \dfrac{\lambda_1^{N-n}+ \delta t p^{n}\lambda_2^{N-n}+\delta t f(I^{n+1}) \lambda_3^{N-n}}{1+\delta t p^{n}+\delta t f(I^{n+1})+\delta t m_s}, \\ \lambda_2^{N-n-1} \gets \dfrac{\lambda_2^{N-n}+\delta t \theta f(I^{n+1})\lambda_3^{N-n}}{1+\delta t \theta f(I^{n+1}) +\delta t m_v}, \\ \lambda_3^{N-n-1} \gets \dfrac{\lambda_3^{N-n}+ \delta t \alpha+ \delta t \varepsilon\lambda_4^{N-n}}{1+ \delta t \varepsilon+ \delta t m_e}, \\ \lambda_4^{N-n-1} \gets \dfrac{\lambda_4^{N-n} + \delta t \beta - \delta t f'(I^{n+1}) S^{n+1} (\lambda_1^{N-n-1}-\lambda_3^{N-n-1}) }{1+ \delta t (m_i+\zeta+\gamma)} \\ \qquad \qquad \quad -\dfrac{\delta t \theta f'(I^{n+1})V^{n+1}(\lambda_2^{N-n-1}- \lambda_3^{N-n-1})}{1+ \delta t (m_i+\zeta+\gamma)} \\ p^{n+1} \gets \max(\min(\dfrac{(\lambda_1^{N-n-1}-\lambda_2^{N-n-1}) S^{n+1}}{a}, p_{\max}), p_{\min}). \end{array} \right. \end{eqnarray*} end |
[1] | CDC, Measles (Rubeola), 2020. Available from: https://www.cdc.gov/measles/index.html. |
[2] |
M. Fakhruddin, D. Suandi, Sumiati, H. Fahlena, N. Nuraini, E. Soewono, Investigation of a measles transmission with vaccination: A case study in Jakarta, Indonesia. Math. Biosci. Eng., 17 (2020), 2998–3018. https://doi.org/10.3934/mbe.2020170. doi: 10.3934/mbe.2020170
![]() |
[3] | WHO, Measles, 2019. Available from: https://www.who.int/en/news-room/fact-sheets/detail/measles. |
[4] | O. Diekmann, J. Heesterbeek, Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis, and Interpretation, Jhon Wiley, 2000. |
[5] |
P. Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Math. Biosci., 180 (2002), 29–48. https://doi.org/10.1016/S0025-5564(02)00108-6. doi: 10.1016/S0025-5564(02)00108-6
![]() |
[6] |
A. A. Alderremy, J. F. Gómez-Aguilar, S. Aly, K. M. Saad, A fuzzy fractional model of coronavirus (COVID-19) and its study with Legendre spectral method, Results Phys., 21 (2021), 103773. https://doi.org/10.1016/j.rinp.2020.103773. doi: 10.1016/j.rinp.2020.103773
![]() |
[7] |
D. Aldila, D. Asrianti, A deterministic model of measles with imperfect vaccination and quarantine intervention, J. Phys.: Conf. Ser., 1218 (2019), 012044. https://doi.org/10.1088/1742-6596/1218/1/012044. doi: 10.1088/1742-6596/1218/1/012044
![]() |
[8] |
M. Sen, S. Alonso-Quesada, Vaccination strategies based on feedback control techniques for a general SEIR-epidemic model, Appl. Math. Comput., 218 (2011), 3888–3904. https://doi.org/10.1016/j.amc.2011.09.036. doi: 10.1016/j.amc.2011.09.036
![]() |
[9] |
M. El Hajji, S. Sayari, Analysis of a fractional-order "SVEIR" epidemic model with a general nonlinear saturated incidence rate in a continuous reactor, Asian Res. J. Math., 12 (2019), 1–17. https://doi.org/10.9734/arjom/2019/v12i430095. doi: 10.9734/arjom/2019/v12i430095
![]() |
[10] |
M. El Hajji, Boundedness and asymptotic stability of nonlinear Volterra integro-differential equations using Lyapunov functional, J. King Saud Univ., Sci., 31 (2019), 1516–1521. https://doi.org/10.1016/j.jksus.2018.11.012. doi: 10.1016/j.jksus.2018.11.012
![]() |
[11] |
M. El Hajji, How can inter-specific interferences explain coexistence or confirm the competitive exclusion principle in a chemostat, Int. J. Biomath., 11 (2018), 1850111. https://doi.org/10.1142/S1793524518501115. doi: 10.1142/S1793524518501115
![]() |
[12] |
A. B. Gumel, C. C. McCluskey, J. Watmough, An sveir model for assessing potential impact of an imperfect anti-sars vaccine, Math. Biosci. Eng., 3 (2006), 485–512. https://doi.org/10.3934/mbe.2006.3.485. doi: 10.3934/mbe.2006.3.485
![]() |
[13] | M. El Hajji, N. Chorfi, M. Jleli, Mathematical modelling and analysis for a three-tiered microbial food web in a chemostat, Electron. J. Differ. Equations, 2017 (2017), 1–13. Available from: http://ejde.math.unt.edu. |
[14] | M. El Hajji, N. Chorfi, M Jleli, Mathematical model for a membrane bioreactor process, Electron. J. Differ. Equations, 2015 (2015), 1–7. Available from: http://ejde.math.txstate.edu. |
[15] |
M. Farman, A. Ahmad, M. U. Saleem, M. O. Ahmad, Analysis and numerical solution of epidemic models by using nonstandard finite difference scheme, Pure Appl. Biol., 9 (2020), 674–682. http://dx.doi.org/10.19045/bspab.2020.90073. doi: 10.19045/bspab.2020.90073
![]() |
[16] |
L. Michel, C. J. Silva, D. F. M. Torres, Model-free based control of a HIV/AIDS prevention model, Math. Biosc. Eng., 19 (2022), 759–774. https://doi.org/10.3934/mbe.2022034. doi: 10.3934/mbe.2022034
![]() |
[17] |
C. J. Silva, G. Cantin, C. Cruz, R. Fonseca-Pinto, R. Passadouro, E. S. dos Santos, et al., Complex network model for COVID-19: Human behavior, pseudo-periodic solutions and multiple epidemic waves, J. Math. Anal. Appli., 2021, 125171. https://doi.org/10.1016/j.jmaa.2021.125171. doi: 10.1016/j.jmaa.2021.125171
![]() |
[18] |
H. W. Hethcote, The mathematics of infectious diseases, SIAM Rev., 42 (2000), 599–653. https://doi.org/10.1137/S0036144500371907. doi: 10.1137/S0036144500371907
![]() |
[19] |
I. A. Moneim, An SEIR model with infectious latent and a periodic vaccination strategy, Math. Modell. Anal., 26 (2021), 236–252. https://doi.org/10.3846/mma.2021.12945. doi: 10.3846/mma.2021.12945
![]() |
[20] |
M. El Hajji, Modelling and optimal control for Chikungunya disease, Theory Biosci., 140 (2021), 27–44. https://doi.org/10.1007/s12064-020-00324-4. doi: 10.1007/s12064-020-00324-4
![]() |
[21] |
M. El Hajji, S. Sayari, A. Zaghdani, Mathematical analysis of an "SIR" epidemic model in a continuous reactor—deterministic and probabilistic approaches, J. Korean Math. Soc., 58 (2021), 45–67. https://doi.org/10.4134/JKMS.j190788. doi: 10.4134/JKMS.j190788
![]() |
[22] |
H. W. Hethcote, Three basic epidemiological models, Biomathematics, 18 (1989). https://doi.org/10.1007/978-3-642-61317-3_5. doi: 10.1007/978-3-642-61317-3_5
![]() |
[23] |
W. O. Kermack, A. G. McKendrick, A contribution to the mathematical theory of epidemics, R. Soc., 115 (1927), 700–721. https://doi.org/10.1098/rspa.1927.0118. doi: 10.1098/rspa.1927.0118
![]() |
[24] |
S. Edward, E. K. Raymond, T. K. Gabriel, F. Nestory, G. M. Godfrey, P. M. Arbogast, A mathematical model for control and elimination of the transmission dynamics of measles, Appl. Comput. Math., 4 (2015), 396–408. http://doi.org/10.11648/j.acm.20150406.12. doi: 10.11648/j.acm.20150406.12
![]() |
[25] |
A. A. Momoh, M. O. Ibrahim, I. J. Uwanta, S. B. Manga, Mathematical model for control of measles epidemiology, I. J. Pure Appl. Math., 87 (2013), 707–718. https://doi.org/10.12732/ijpam.v87i5.4. doi: 10.12732/ijpam.v87i5.4
![]() |
[26] |
H. Wei, Y. Jiang, X. Song, G. H. Su, S. Z. Qiu, Global attractivity and permanence of a SVEIR epidemic model with pulse vaccination and time delay, J. Comput. Appl. Math., 229 (2009), 302–312. https://doi.org/10.1016/j.cam.2008.10.046. doi: 10.1016/j.cam.2008.10.046
![]() |
[27] | P. Adda, L. N. Nkague, G. Sallet, L. Castelli, A SVEIR model with imperfect vaccine, in CMPD 3 Conference on Computational and Mathematical Population Dynamics, 2010. https://hal.inria.fr/hal-00764764. |
[28] |
L. N. Nkague, J. M. Ntaganda, H. Abboubakar, J. C. Kamgang, L. Castelli, Global stability of a SVEIR epidemic model: Application to poliomyelitis transmission dynamics, Open J. Modell. Simul., 5 (2017), 98–112. https://doi.org/10.4236/ojmsi.2017.51008. doi: 10.4236/ojmsi.2017.51008
![]() |
[29] | J. P. LaSalle, The Stability of Dynamical Systems, SIAM, 25 (1976). https://doi.org/10.1137/1.9781611970432. |
[30] |
M. El Hajji, A. Zaghdani, S. Sayari, Mathematical analysis and optimal control for Chikungunya virus with two routes of infection with nonlinear incidence rate, Int. J. Biomath., (2021), 2150088. https://doi.org/10.1142/S1793524521500881. doi: 10.1142/S1793524521500881
![]() |
[31] |
N. Chitnis, J. M. Hyman, J. M. Cushing, Determining important parameters in the spread of malaria through the sensitivity analysis of a mathematical model, Bull. Math. Biol., 70 (2008), 1272–1296. https://doi.org/10.1007/s11538-008-9299-0. doi: 10.1007/s11538-008-9299-0
![]() |
[32] |
C. J. Silva, D. F. M. Torres, Optimal control for a tuberculosis model with reinfection and post-exposure interventions, Math. Biosci., 244 (2013), 154–164. https://doi.org/10.1016/j.mbs.2013.05.005. doi: 10.1016/j.mbs.2013.05.005
![]() |
[33] |
H. S. Rodrigues, M. T. T. Monteiro, D. F. M. Torres, Sensitivity analysis in a dengue epidemiological model, Conf. Pap. Sci., 2013. https://doi.org/10.1155/2013/721406. doi: 10.1155/2013/721406
![]() |
[34] | W. Fleming, R. Rishel, Deterministic and Stochastic Optimal Control, Springer Verlag, New York, 1975. https://doi.org/10.1007/978-1-4612-6380-7. |
[35] | S. Lenhart, J. T. Workman, Optimal Control Applied to Biological Models, Chapman and Hall, 2007. https://doi.org/10.1201/9781420011418. |
[36] | L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, E. F. Mishchenko, K. N. Trirogoff, L. W. Neustadt, The Mathematical Theory of Optimal Processes, 2000. https://doi.org/10.1201/9780203749319. |
[37] |
A. B. Gumel, P. N. Shivakumar, B. M. Sahai, A mathematical model for the dynamics of HIV-1 during the typical course of infection, Nonlinear Anal., Theory Methods Appl., 47 (2001), 1773–1783. https://doi.org/10.1016/S0362-546X(01)00309-1. doi: 10.1016/S0362-546X(01)00309-1
![]() |
1. | Mourad Chamekh, Mohamed Ali Latrach, Fadel Jday, Multi-step semi-analytical solutions for a chikungunya virus system, 2023, 2731-6734, 10.1007/s43994-023-00027-8 | |
2. | Muhammad Farman, Aamir Shehzad, Ali Akgül, Dumitru Baleanu, Manuel De la Sen, Modelling and Analysis of a Measles Epidemic Model with the Constant Proportional Caputo Operator, 2023, 15, 2073-8994, 468, 10.3390/sym15020468 | |
3. | Olumuyiwa James Peter, Sania Qureshi, Mayowa M. Ojo, Ratchada Viriyapong, Amanullah Soomro, Mathematical dynamics of measles transmission with real data from Pakistan, 2022, 2363-6203, 10.1007/s40808-022-01564-7 | |
4. | Salah Alsahafi, Stephen Woodcock, Exploring HIV Dynamics and an Optimal Control Strategy, 2022, 10, 2227-7390, 749, 10.3390/math10050749 | |
5. | Zakaria Yaagoub, Karam Allali, Global Stability of Multi-Strain SEIR Epidemic Model with Vaccination Strategy, 2023, 28, 2297-8747, 9, 10.3390/mca28010009 | |
6. | Ahmed Alshehri, Miled El Hajji, Mathematical study for Zika virus transmission with general incidence rate, 2022, 7, 2473-6988, 7117, 10.3934/math.2022397 | |
7. | Nikhila Yaladanda, Rajasekhar Mopuri, Hari Prasad Vavilala, Srinivasa Rao Mutheneni, Modelling the impact of perfect and imperfect vaccination strategy against SARS CoV-2 by assuming varied vaccine efficacy over India, 2022, 15, 22133984, 101052, 10.1016/j.cegh.2022.101052 | |
8. | Mohammed H. Alharbi, Global investigation for an "SIS" model for COVID-19 epidemic with asymptomatic infection, 2023, 20, 1551-0018, 5298, 10.3934/mbe.2023245 | |
9. | Amer Hassan Albargi, Miled El Hajji, Mathematical analysis of a two-tiered microbial food-web model for the anaerobic digestion process, 2023, 20, 1551-0018, 6591, 10.3934/mbe.2023283 | |
10. | Yantao Luo, Yunqiang Yuan, Zhidong Teng, ANALYSIS OF A DEGENERATED DIFFUSION SVEQIRV EPIDEMIC MODEL WITH GENERAL INCIDENCE IN A SPACE HETEROGENEOUS ENVIRONMENT, 2024, 14, 2156-907X, 2704, 10.11948/20230346 | |
11. | Rukhsar Ikram, Amir Khan, Mostafa Zahri, Aeshah A. Raezah, The impact of Lévy noise on a stochastic measles model, 2023, 138, 2190-5444, 10.1140/epjp/s13360-023-04625-7 | |
12. | Maoxing Liu, Yuhang Li, Dynamics analysis of an SVEIR epidemic model in a patchy environment, 2023, 20, 1551-0018, 16962, 10.3934/mbe.2023756 | |
13. | Dilber Uzun Ozsahin, Najeeb Alam Khan, Araib Aqeel, Hijaz Ahmad, Maged F. Alotaibi, Muhammad Ayaz, Junyuan Yang, Mathematical modeling and dynamics of immunological exhaustion caused by measles transmissibility interaction with HIV host, 2024, 19, 1932-6203, e0297476, 10.1371/journal.pone.0297476 | |
14. | Liang’an 良安 Huo 霍, Yue 跃 Yu 于, Impact of individual behavior adoption heterogeneity on epidemic transmission in multiplex networks, 2023, 32, 1674-1056, 108703, 10.1088/1674-1056/acea65 | |
15. | Jiandong Nie, Qiaoling Chen, Zhidong Teng, Yihan Zhang, Ramziya Rifhat, Dynamics of a Stochastic Measles Model with General Incidence Rate and Black–Karasinski Process, 2024, 47, 0126-6705, 10.1007/s40840-024-01771-8 | |
16. | Miled El Hajji, Periodic solutions for chikungunya virus dynamics in a seasonal environment with a general incidence rate, 2023, 8, 2473-6988, 24888, 10.3934/math.20231269 | |
17. | Miled El Hajji, Dalal M. Alshaikh, Nada A. Almuallem, Periodic Behaviour of an Epidemic in a Seasonal Environment with Vaccination, 2023, 11, 2227-7390, 2350, 10.3390/math11102350 | |
18. | Young Rock Kim, Youngho Min, Joy Nana Okogun-Odompley, Rehana Naz, A mathematical model of COVID-19 with multiple variants of the virus under optimal control in Ghana, 2024, 19, 1932-6203, e0303791, 10.1371/journal.pone.0303791 | |
19. | Marlon Nunes Gonzaga, Marcelo Martins de Oliveira, Allbens Picardi Faria Atman, Role of Vaccination Strategies to Host-Pathogen Dynamics in Social Interactions, 2024, 26, 1099-4300, 739, 10.3390/e26090739 | |
20. | Amar Nath Chatterjee, Santosh Kumar Sharma, Fahad Al Basir, A Compartmental Approach to Modeling the Measles Disease: A Fractional Order Optimal Control Model, 2024, 8, 2504-3110, 446, 10.3390/fractalfract8080446 | |
21. | Amer Hassan Albargi, Miled El Hajji, Bacterial Competition in the Presence of a Virus in a Chemostat, 2023, 11, 2227-7390, 3530, 10.3390/math11163530 | |
22. | W. Ahmad, A. I. K. Butt, M. Rafiq, Z. Asif, T. Ismaeel, N. Ahmad, Modeling, analyzing and simulating the Measles transmission dynamics through efficient computational optimal control technique, 2024, 139, 2190-5444, 10.1140/epjp/s13360-024-05368-9 | |
23. | Zakaria Yaagoub, Fractional two-strain SVLIR epidemic model with vaccination and quarantine strategies, 2025, 13, 2195-268X, 10.1007/s40435-024-01561-x | |
24. | Tianqi Song, Yishi Wang, Chuncheng Wang, Qi An, Modeling and analyzing COVID-19 infections in South Africa, 2025, 18, 1793-5245, 10.1142/S1793524523501036 | |
25. | Nassira Madani, Zakia Hammouch, EL-Houssine Azroul, New model of HIV/AIDS dynamics based on Caputo–Fabrizio derivative order: Optimal strategies to control the spread, 2025, 18777503, 102612, 10.1016/j.jocs.2025.102612 |
Variable | Description | Parameter | Description |
S(t) | Susceptible individuals | f | Saturated incidence rate |
V(t) | Vaccinated individuals | \Lambda | Susceptible recruitment rate |
E(t) | Exposed individuals | \varepsilon | Infection transmission rate |
I(t) | Infected individuals | \gamma | Recovery rate |
R(t) | Recovered individuals | \theta\in [0, 1] | Effectiveness rate of the vaccination |
Parameter | Description | Parameter | Description |
m_v | Mortality rate of vaccinated individuals | m_s | mortality rate of susceptible individuals |
m_e | Mortality rate of exposed individuals | p | vaccination rate |
m_i+\zeta | Total mortality rate of infected individuals | \dfrac{1}{\varepsilon} | Average latency time spent in compartment E |
m_r | Mortality rate of recovered individuals | \dfrac{1}{\gamma} | Average duration elapsed in compartment I |
Parameter | Sensitivity index of \mathcal{R} | Sign |
\Lambda | Y^{ \mathcal{R}}_{\Lambda}=\dfrac{\partial \mathcal{R}}{\partial \Lambda} \times \dfrac{\Lambda}{ \mathcal{R}}=1 | +ve |
\varepsilon | Y^{ \mathcal{R}}_{\varepsilon}=\dfrac{\partial \mathcal{R}}{\partial \varepsilon} \times \dfrac{\varepsilon}{ \mathcal{R}}=\dfrac{m_e}{m_e+\varepsilon} | +ve |
\theta | Y^{ \mathcal{R}}_{\theta}=\dfrac{\partial \mathcal{R}}{\partial \theta} \times \dfrac{\theta}{ \mathcal{R}}=\dfrac{\theta p}{m_v+\theta p} | +ve |
p | Y^{ \mathcal{R}}_{p}=\dfrac{\partial \mathcal{R}}{\partial p} \times \dfrac{p}{ \mathcal{R}}=\dfrac{(\theta m_s- m_v)p}{(m_s+ p)(m_v+\theta p)} | |
m_v | Y^{ \mathcal{R}}_{m_v}=\dfrac{\partial \mathcal{R}}{\partial m_v} \times \dfrac{m_v}{ \mathcal{R}}=\dfrac{-\theta p}{m_v+\theta p} | -ve |
m_s | Y^{ \mathcal{R}}_{m_s}=\dfrac{\partial \mathcal{R}}{\partial m_s} \times \dfrac{m_s}{ \mathcal{R}}=-\dfrac{m_s}{m_s+p} | -ve |
m_e | Y^{ \mathcal{R}}_{m_e}=\dfrac{\partial \mathcal{R}}{\partial m_e} \times \dfrac{m_e}{ \mathcal{R}}=-\dfrac{m_e}{m_e+\varepsilon} | -ve |
m_i | Y^{ \mathcal{R}}_{m_i}=\dfrac{\partial \mathcal{R}}{\partial m_i} \times \dfrac{m_i}{ \mathcal{R}}=-\dfrac{m_i}{m_i+\zeta+\gamma} | -ve |
\zeta | Y^{ \mathcal{R}}_{\zeta}=\dfrac{\partial \mathcal{R}}{\partial \zeta} \times \dfrac{\zeta}{ \mathcal{R}}=-\dfrac{\zeta}{m_i+\zeta+\gamma} | -ve |
\gamma | Y^{ \mathcal{R}}_{\gamma}=\dfrac{\partial \mathcal{R}}{\partial \gamma} \times \dfrac{\gamma}{ \mathcal{R}}=-\dfrac{\gamma}{m_i+\zeta+\gamma} | -ve |
Parameter | \bar f | k | m_s | m_v | m_e | m_i | \zeta | m_r | p | \theta |
Value | 5 | 3 | 0.2 | 0.1 | 0.3 | 0.25 | 0.15 | 0.2 | 1 | 0.5 |
Algorithm 1: Optimal control resolution |
S^0 \gets S_0, V^0 \gets V_0, E^0 \gets E_0 , I^0 \gets I_0 , \lambda_1^N \gets 0, \lambda_2^N \gets 0, \lambda_3^N \gets 0, \lambda_4^N \gets 0 , for n = 0 to N-1 do \begin{eqnarray*} \left[ \begin{array}{llll} S^{n+1} \gets \dfrac{S^n+\delta t \Lambda}{1+\delta t (m_s+p^n)+ \delta tf(I^n)}, \\ V^{n+1} \gets \dfrac{V^n+\delta t p^nS^{n+1}}{1+\delta t m_v+ \delta t\theta f(I^n)}, \\ E^{n+1} \gets \dfrac{E^n+\delta t f(I^n)(S^{n+1}+\theta V^{n+1})}{1+ \delta t (m_e+\varepsilon)}, \\ I^{n+1} \gets \dfrac{I^n+\delta t \varepsilon E^{n+1}}{1+\delta t (m_i+\zeta+\gamma)}, \\ \lambda_1^{N-n-1} \gets \dfrac{\lambda_1^{N-n}+ \delta t p^{n}\lambda_2^{N-n}+\delta t f(I^{n+1}) \lambda_3^{N-n}}{1+\delta t p^{n}+\delta t f(I^{n+1})+\delta t m_s}, \\ \lambda_2^{N-n-1} \gets \dfrac{\lambda_2^{N-n}+\delta t \theta f(I^{n+1})\lambda_3^{N-n}}{1+\delta t \theta f(I^{n+1}) +\delta t m_v}, \\ \lambda_3^{N-n-1} \gets \dfrac{\lambda_3^{N-n}+ \delta t \alpha+ \delta t \varepsilon\lambda_4^{N-n}}{1+ \delta t \varepsilon+ \delta t m_e}, \\ \lambda_4^{N-n-1} \gets \dfrac{\lambda_4^{N-n} + \delta t \beta - \delta t f'(I^{n+1}) S^{n+1} (\lambda_1^{N-n-1}-\lambda_3^{N-n-1}) }{1+ \delta t (m_i+\zeta+\gamma)} \\ \qquad \qquad \quad -\dfrac{\delta t \theta f'(I^{n+1})V^{n+1}(\lambda_2^{N-n-1}- \lambda_3^{N-n-1})}{1+ \delta t (m_i+\zeta+\gamma)} \\ p^{n+1} \gets \max(\min(\dfrac{(\lambda_1^{N-n-1}-\lambda_2^{N-n-1}) S^{n+1}}{a}, p_{\max}), p_{\min}). \end{array} \right. \end{eqnarray*} end |
Variable | Description | Parameter | Description |
S(t) | Susceptible individuals | f | Saturated incidence rate |
V(t) | Vaccinated individuals | \Lambda | Susceptible recruitment rate |
E(t) | Exposed individuals | \varepsilon | Infection transmission rate |
I(t) | Infected individuals | \gamma | Recovery rate |
R(t) | Recovered individuals | \theta\in [0, 1] | Effectiveness rate of the vaccination |
Parameter | Description | Parameter | Description |
m_v | Mortality rate of vaccinated individuals | m_s | mortality rate of susceptible individuals |
m_e | Mortality rate of exposed individuals | p | vaccination rate |
m_i+\zeta | Total mortality rate of infected individuals | \dfrac{1}{\varepsilon} | Average latency time spent in compartment E |
m_r | Mortality rate of recovered individuals | \dfrac{1}{\gamma} | Average duration elapsed in compartment I |
Parameter | Sensitivity index of \mathcal{R} | Sign |
\Lambda | Y^{ \mathcal{R}}_{\Lambda}=\dfrac{\partial \mathcal{R}}{\partial \Lambda} \times \dfrac{\Lambda}{ \mathcal{R}}=1 | +ve |
\varepsilon | Y^{ \mathcal{R}}_{\varepsilon}=\dfrac{\partial \mathcal{R}}{\partial \varepsilon} \times \dfrac{\varepsilon}{ \mathcal{R}}=\dfrac{m_e}{m_e+\varepsilon} | +ve |
\theta | Y^{ \mathcal{R}}_{\theta}=\dfrac{\partial \mathcal{R}}{\partial \theta} \times \dfrac{\theta}{ \mathcal{R}}=\dfrac{\theta p}{m_v+\theta p} | +ve |
p | Y^{ \mathcal{R}}_{p}=\dfrac{\partial \mathcal{R}}{\partial p} \times \dfrac{p}{ \mathcal{R}}=\dfrac{(\theta m_s- m_v)p}{(m_s+ p)(m_v+\theta p)} | |
m_v | Y^{ \mathcal{R}}_{m_v}=\dfrac{\partial \mathcal{R}}{\partial m_v} \times \dfrac{m_v}{ \mathcal{R}}=\dfrac{-\theta p}{m_v+\theta p} | -ve |
m_s | Y^{ \mathcal{R}}_{m_s}=\dfrac{\partial \mathcal{R}}{\partial m_s} \times \dfrac{m_s}{ \mathcal{R}}=-\dfrac{m_s}{m_s+p} | -ve |
m_e | Y^{ \mathcal{R}}_{m_e}=\dfrac{\partial \mathcal{R}}{\partial m_e} \times \dfrac{m_e}{ \mathcal{R}}=-\dfrac{m_e}{m_e+\varepsilon} | -ve |
m_i | Y^{ \mathcal{R}}_{m_i}=\dfrac{\partial \mathcal{R}}{\partial m_i} \times \dfrac{m_i}{ \mathcal{R}}=-\dfrac{m_i}{m_i+\zeta+\gamma} | -ve |
\zeta | Y^{ \mathcal{R}}_{\zeta}=\dfrac{\partial \mathcal{R}}{\partial \zeta} \times \dfrac{\zeta}{ \mathcal{R}}=-\dfrac{\zeta}{m_i+\zeta+\gamma} | -ve |
\gamma | Y^{ \mathcal{R}}_{\gamma}=\dfrac{\partial \mathcal{R}}{\partial \gamma} \times \dfrac{\gamma}{ \mathcal{R}}=-\dfrac{\gamma}{m_i+\zeta+\gamma} | -ve |
Parameter | \bar f | k | m_s | m_v | m_e | m_i | \zeta | m_r | p | \theta |
Value | 5 | 3 | 0.2 | 0.1 | 0.3 | 0.25 | 0.15 | 0.2 | 1 | 0.5 |
Algorithm 1: Optimal control resolution |
S^0 \gets S_0, V^0 \gets V_0, E^0 \gets E_0 , I^0 \gets I_0 , \lambda_1^N \gets 0, \lambda_2^N \gets 0, \lambda_3^N \gets 0, \lambda_4^N \gets 0 , for n = 0 to N-1 do \begin{eqnarray*} \left[ \begin{array}{llll} S^{n+1} \gets \dfrac{S^n+\delta t \Lambda}{1+\delta t (m_s+p^n)+ \delta tf(I^n)}, \\ V^{n+1} \gets \dfrac{V^n+\delta t p^nS^{n+1}}{1+\delta t m_v+ \delta t\theta f(I^n)}, \\ E^{n+1} \gets \dfrac{E^n+\delta t f(I^n)(S^{n+1}+\theta V^{n+1})}{1+ \delta t (m_e+\varepsilon)}, \\ I^{n+1} \gets \dfrac{I^n+\delta t \varepsilon E^{n+1}}{1+\delta t (m_i+\zeta+\gamma)}, \\ \lambda_1^{N-n-1} \gets \dfrac{\lambda_1^{N-n}+ \delta t p^{n}\lambda_2^{N-n}+\delta t f(I^{n+1}) \lambda_3^{N-n}}{1+\delta t p^{n}+\delta t f(I^{n+1})+\delta t m_s}, \\ \lambda_2^{N-n-1} \gets \dfrac{\lambda_2^{N-n}+\delta t \theta f(I^{n+1})\lambda_3^{N-n}}{1+\delta t \theta f(I^{n+1}) +\delta t m_v}, \\ \lambda_3^{N-n-1} \gets \dfrac{\lambda_3^{N-n}+ \delta t \alpha+ \delta t \varepsilon\lambda_4^{N-n}}{1+ \delta t \varepsilon+ \delta t m_e}, \\ \lambda_4^{N-n-1} \gets \dfrac{\lambda_4^{N-n} + \delta t \beta - \delta t f'(I^{n+1}) S^{n+1} (\lambda_1^{N-n-1}-\lambda_3^{N-n-1}) }{1+ \delta t (m_i+\zeta+\gamma)} \\ \qquad \qquad \quad -\dfrac{\delta t \theta f'(I^{n+1})V^{n+1}(\lambda_2^{N-n-1}- \lambda_3^{N-n-1})}{1+ \delta t (m_i+\zeta+\gamma)} \\ p^{n+1} \gets \max(\min(\dfrac{(\lambda_1^{N-n-1}-\lambda_2^{N-n-1}) S^{n+1}}{a}, p_{\max}), p_{\min}). \end{array} \right. \end{eqnarray*} end |