A mathematical investigation of an "SVEIR" epidemic model for the measles transmission

1.   Introduction

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 $\mathcal{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 $\mathcal{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.

2.   Mathematical model and results

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.

● $\Lambda$ is the recruitment rate of susceptible individuals.

● $p$ describe the rate of vaccination.

● $\theta\in [0, 1]$ is the effectiveness rate of the vaccination.

● Each compartment has a specific natural mortality rate. $m_s$, $m_v$, $m_e$, $m_i$ and $m_r$ are the mortality rates of susceptible, vaccinated, exposed, infected and recovered individuals, respectively.

● Measles can cause death (fatal) with a mortality rate of $\zeta$.

● $f$ represents the saturated incidence rate.

● $\varepsilon$ the rate at which an exposed individuals become infectious. $1/\varepsilon$ is the average latency time spent in compartment $E$.

● $\gamma$ the rate at which infectious agents recover their health. $1/\gamma$ 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 $\varepsilon$.

● The recovery rate of individuals infected with measles is $\gamma$.

● 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:

with positive initial condition $(S(0), V(0), E(0), I(0), R(0))\in \mathbb{R}_+^5\;$.

Assumption 1. The function $f$ is non-negative $C^1(\mathbb{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) \leq f'(0) I$ for all $I > 0$.

Proof. For $I \in \mathbb{R}_+$, let $h_1(I) = f(I)-If'(I)$. Since $f$ is an increasing and concave function then $f'(I) \geq 0$ and $f''(I) \leq 0$. Therefore $h_1'(I) = -If''(I) \geq 0$ and $h_1(I) \geq h_1(0) = 0$ or also $f(I)\geq If'(I)$. Similarly, let $h_2(I) = f(I)-If'(0)$ then $h_2'(I) = f'(I)-f'(0) \leq 0$ since $f$ is concave. Then $h_2(I) \leq h_2(0) = 0$ then $f(I)\leq If'(0)$.

Remark 1. The Monod (or Holling's type Ⅱ) functions satisfy the assumption 1 and then can express the transmission rate.

$\bar 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(m_s, m_v, m_e, m_i+\zeta, m_r)$, then

Proposition 1. The compact set $\Omega_1 = \{(S, V, E, I, R)\in \mathbb{R}_+^5\; /\; S+V+E+I+R\leq \dfrac{\Lambda}{m} \}$ is positively invariant for system (2.1).

Proof. Since $\dot S_{|S = 0} = \Lambda > 0$, $\dot V_{|V = 0} = pS > 0$, $\dot E_{|E = 0} = f(I)(S+\theta V) > 0$, $\dot I_{|I = 0} = \varepsilon E > 0$ and $\dot R_{|R = 0} = \gamma I > 0$, then the solution of system (2.1) is non-negative.

Consider a new variable $T = S+V+E+I+R-\dfrac{\Lambda}{m}$, then by adding all equations of (2.1), $T$ satisfies:

Hence

Since all variables are non-negative then all variables are bounded and thus $\Omega_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 $\mathcal{R}$. Recall that $\mathcal{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 $\mathcal{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 $\mathcal{R}$. Let $S_0 = \dfrac{\Lambda}{m_s+p}, V_0 = \dfrac{p \Lambda}{m_v (m_s+p)}, E_0 = 0, \bar I_0 = 0$ and $R_0 = 0$. The Jacobian of the sub-system $(E, I)$ of the original system (2.1) at the disease-free equilibrium point $\mathcal{E}_0 = (S_0, V_0, E_0, I_0, R_0)$ is given by

where $F = \left(\begin{array}{ccccccc} 0 &f'(0)(S_0+\theta V_0)\\ 0 & 0 \end{array} \right)$ and $V = \left(\begin{array}{ccccccc} (m_e+\varepsilon) & 0\\ -\varepsilon & (m_i+\zeta +\gamma) \end{array} \right)$.

The determinant of $V$ is given by det$(V) = (m_e+\varepsilon)(m_i+\zeta +\gamma) > 0$ and thus the inverse matrix of $V$ is

and the next-generation matrix is given by

Then, the basic reproduction number of system (2.1) is calculated as the spectral radius of the matrix $FV^{-1}$:

Note that in absence of vaccination ($p = 0$), the system (2.1) becomes the standard SEIR model with

Proposition 2. System (2.1) admits a trivial (disease-free) equilibrium point $\mathcal{E}_0 = (S_0, V_0, E_0, I_0, R_0)$ and if $\mathcal{R} > 1$, system (2.1) admits a unique non-trivial (endemic) equilibrium point $\mathcal{E}_1 = (S_1, V_1, E_1, I_1, R_1)$ with $S_1, V_1, E_1, I_1, R_1 > 0$.

Proof. An equilibrium point for the system (2.1) satisfies

which reduces to

From the fourth equation of (2.6) one deduces that

Since all parameters are non-negative then either $I = 0$ or

● If $I = 0$ then $R = 0, E = 0, S = \dfrac{\Lambda}{m_s+p} = \bar S$ and $V = \dfrac{pS}{m_v} = \dfrac{p \Lambda}{m_v(m_s+p)} = \bar V$. This equilibrium known as the disease-free equilibrium denoted here by $\mathcal{E}_0 = (S_0, V_0, E_0, I_0, R_0)$.

● For $I\not = 0$, define the continuous function $g$

The derivative of $g$ is given by

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

and

Since $\mathcal{R} > 1$, $\lim_{I \mapsto 0}g(I) < 0$ and $g(\dfrac{\Lambda}{m}) > 0$, the equation $g(I) = 0$ admits a unique solution $I_1$ in $(0, \dfrac{\Lambda}{m})$ and then the uniqueness of the endemic equilibrium $\mathcal{E}_1 = (S_1, V_1, E_1, I_1, R_1)$.

In the section below, we discussed the local stability of equilibria with respect to the value of the basic reproduction number $\mathcal{R}$.

3.   Local stability analysis of the DFE

Theorem 1. If $\mathcal{R} < 1$, then the disease-free equilibrium (DFE) point $\mathcal{E}_0$ is locally asymptotically stable and it is unstable if $\mathcal{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:

It's value, evaluated at $\mathcal{E}_0$, is given by:

which admits five distinguish eigenvalues given by $\lambda_1 = -(m_s+p) < 0, \lambda_2 = -m_v < 0, \lambda_3 = -m_r < 0$ and two other eigenvalues corresponding to the sub-matrix

The corresponding characteristic polynomial is calculated as follows

Only if $\mathcal{R} < 1$, then the roots of $P(.)$ have negative real parts. Then $\mathcal{E}_0$ is locally asymptotically stable only if $\mathcal{R} < 1$. $\mathcal{E}_0$ is unstable once $\mathcal{R} > 1$.

In the next section, the global stability behavior of the equilibria will be discussed.

4.   Global stability analysis

Corollary 1. The set $\Omega_2 = \{(S, V, E, I, R)\in \mathbb{R}_+^5\; /\; S+V+E+I+R\leq \dfrac{\Lambda}{m}; S\leq S_0, V\leq V_0\}$ is positively invariant for system (2.1).

Proof. It is already proved that $\Omega_1$ is invariant for system (2.1). Note that $\dot S(t) < 0$ for $S(t) > S_0$ therefore $\liminf S(t)\leq S_0$. Let $\xi > 0$ be a constant and $S(0)$ be an initial condition. Then $\exists \; T\geq 0; \; S(t)\leq S_0+\xi\; \forall t\geq T$. Therefore

Consider $\bar V = V -\dfrac{p(S_0+\xi)}{m_v}$ which satisfies

Then

Therefore $\liminf \bar V(t)\leq 0$ and $\liminf V(t)\leq \dfrac{p(S_0+\xi)}{m_v}$ for all $\xi > 0$ thus $\liminf V(t)\leq \dfrac{p S_0}{m_v} = V_0$. This completes the proof.

Theorem 2. $\mathcal{E}_0$ is globally asymptotically stable if $\mathcal{R} \leq 1$ however it is unstable if $\mathcal{R} > 1$.

Proof. Let the function $\mathrm{F}_0$ given by:

Since $f(I) \leq f'(0) I$, the time-derivative of $\mathrm{F}_0$ is given by:

If $\mathcal{R}\leq1,$ then ${\dot {\mathrm{F}}}_0\leq 0$ for all $S, V, E, I, R > 0.$ Let $W_0 = \{(S, V, E, I, R):{\dot {\mathrm{F}}}_0 = 0\} = \{ \mathcal{E}_0\}$. Then by LaSalle's invariance principle ^[29] $\mathcal{E}_0$ is globally asymptotically stable once $\mathcal{R}\leq 1$ (see ^[9,20,21,30] for some applications). Then the solution of system (2.1) converges to $\mathcal{E}_0$ as $t \to +\infty$.

Hereafter, we study the global stability analysis for the endemic equilibrium $\mathcal{E}_1$.

Theorem 3. If $\mathcal{R} > 1$, then the endemic equilibrium $\mathcal{E}_1$ is globally asymptotically stable.

Proof. Let the function $\mathrm{F}_1$, given by

The function $\mathrm{F}_1$ admits its minimum value $\mathrm{F}_{1min} = S_1+V_1+E_1+\dfrac{m_e+\varepsilon}{\varepsilon}I_1$ once $S = S_1, V = V_1, E = E_1, I = I_1$. The time-derivative of $\mathrm{F}_1$, along solutions of system (2.1) is given by

Since

then $(2-\dfrac{S}{S_1}-\dfrac{S_1}{S})\leq 0, (3-\dfrac{S_1}{S}-\dfrac{V}{V_1}-\dfrac{V_1}{V}\dfrac{S}{S_1}) \leq 0,$ $(3-\dfrac{S_1}{S}-\dfrac{I}{I_1}\dfrac{S}{S_1}\dfrac{E_1}{E}-\dfrac{E}{E_1}\dfrac{I_1}{I})\leq 0,$ $(4-\dfrac{V_1}{V}\dfrac{S}{S_1}-\dfrac{I}{I_1}\dfrac{V}{V_1}\dfrac{E_1}{E}-\dfrac{E}{E_1}\dfrac{I_1}{I}-\dfrac{S_1}{S})\leq 0$. Therefore ${\dot {\mathrm{F}}}_2\leq 0$. Thank's to Lyapunov theorem, $\mathcal{E}_1$ is stable.

Now, we have to show the asymptotic stability of $\mathcal{E}_1$ using the Lasalle invariance principle.

Define $A = (2-\dfrac{S}{S_1}-\dfrac{S_1}{S}), B = (3-\dfrac{S_1}{S}-\dfrac{V}{V_1}-\dfrac{V_1}{V}\dfrac{S}{S_1}),$ $C = (3-\dfrac{S_1}{S}-\dfrac{I}{I_1}\dfrac{S}{S_1}\dfrac{E_1}{E}-\dfrac{E}{E_1}\dfrac{I_1}{I})$,

$D = (4-\dfrac{V_1}{V}\dfrac{S}{S_1}-\dfrac{I}{I_1}\dfrac{V}{V_1}\dfrac{E_1}{E}-\dfrac{E}{E_1}\dfrac{I_1}{I}-\dfrac{S_1}{S})$. Then ${\dot {\mathrm{F}}}_1(S, V, E, I, R) = 0$ if and only if $A = B = C = D = 0$.

$A = 0$ means that $S = S_1$ and $B = 0$ means that $V = V_1$. If moreover $C = 0$ then $\dfrac{E}{E_1} = \dfrac{I}{I_1}$. To conclude

Let $a = \dfrac{E}{E_1} = \dfrac{I}{I_1}$, then $E = a E_1$ and $I = a I_1$.

The endemic equilibrium satisfies

Then

and thus $I = I_1$ and $E = E_1$. Finally

Therefore $\{(S, V, E, I, R)\, | {\dot {\mathrm{F}}}_1 = 0\} = \{ \mathcal{E}_1\}$. Then by LaSalle's invariance principle ^[29] $\mathcal{E}_1$ is globally asymptotically stable once $\mathcal{R} > 1$ (see ^[9,20,21,30] for some applications).

5.   Sensitivity analysis

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 $\mathcal{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 $\mathcal{R}$ that depends differentiably on a parameter $\rho$ is defined by

For instance, $Y^{ \mathcal{R}}_{\rho} = 1$ implies an increase (decrease) of $\rho$ by $y\%$ increases (decreases) $\mathcal{R}$ by the same percentage. On the other hand, $Y^{ \mathcal{R}}_{\rho} = -1$ indicates an increase (decrease) of $\rho$ by $y\%$ decreases (increases) $\mathcal{R}$ by $y\%$. In this case, $\rho$ is called a highly sensitive parameter.

Proposition 3. The explicit expression of $\mathcal{R}$ is given by

The basic reproduction number $\mathcal{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 $\Lambda, \varepsilon$ and $\theta$ are positive and hence an increase of the values of parameters $\Lambda, \varepsilon$ and $\theta$ leads to an increase of the value of $\mathcal{R}$. Note that the dependence of $\mathcal{R}$ on $p$ depend on the sign of $(\theta m_s- m_v)$, and the dependence of $\mathcal{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 $\mathcal{R}$.

6.   Optimal control problem

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

where $p = p(t)$ is a function that will be considered as control.

Note that the compact set $\Omega_1 = \{(S, V, E, I)\in \mathbb{R}_+^4\; /\; S+V+E+I\leq \dfrac{\Lambda}{m} \}$ is positively invariant for system (6.1).

Assume moreover that $f$ is globally Lipschitz with an upper bound $\bar f = \sup_{I > 0} f(I)$ and a Lipschitz constant $L_f$. Consider the admissible space

where $p_{\min}$ and $p_{\max}$ are the bounds of the control.

The optimal control problem considered here is

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 $\alpha$ and $\beta$ denote the weight constants of the exposed and infective human population.

For $\varphi = \left(S, V, E, I\right)^t$, the system (6.1) takes a new form

where $A = \left(\begin{array}{cccccc} -m_s-p & 0 & 0 & 0\\ p & -m_v & 0 & 0\\ 0 & 0 & -m_e-\varepsilon & 0 \\ 0 & 0 & \varepsilon & -m_i-\zeta-\gamma \end{array} \right)$ and $G(\varphi) = \left(\begin{array}{ccc} \Lambda-f(I) S\\ -\theta f(I) V \\ (S+\theta V)f(I) \\ 0 \end{array} \right)$.

Proposition 4. The continuous function $F$ is uniformly Lipschitz.

Proof. The continuous function $G$ is uniformly Lipschitz since

where $M = 2 \max(\dfrac{\Lambda}{m} L_f, \theta \dfrac{\Lambda}{m} L_f, \bar f)$. Since

where $\|.\|$ is the matrix norm, then

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)

For a given optimal control $p^*$, there exist adjoint functions $\lambda_1, \lambda_2, \lambda_3$ and $\lambda_4$ corresponding to the states $S, V, E$ and $I$ such that:

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,

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

7.   Numerical simulations

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.

7.1.   Direct problem

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$.

7.2.   Optimal control problems (6.1) and (6.2)

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).

8.   Conclusions

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.

Acknowledgments

The author would like to thank the editors and the anonymous reviewers whose invaluable comments and suggestions have greatly improved this manuscript.

Conflict of interest

The authors declare no conflict of interest.

A.   Appendix

A.1.   Appropriated numerical scheme

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:

Hence, the following algorithm was applied.