Optimal Impulse Vaccination Approach for an SIR Control Model with Short-Term Immunity

Abstract

Vaccines are not administered on a continuous basis, but injections are practically introduced at discrete times often separated by an important number of time units, and this differs depending on the nature of the epidemic and its associated vaccine. In addition, especially when it comes to vaccination, most optimization approaches in the literature and those that have been subject to epidemic models have focused on treating problems that led to continuous vaccination schedules but their applicability remains debatable. In search of a more realistic methodology to resolve this issue, a control modeling design, where the control can be characterized analytically and then optimized, can definitely help to find an optimal regimen of pulsed vaccinations. Therefore, we propose a susceptible-infected-removed (SIR) hybrid epidemic model with impulse vaccination control and a compartment that represents the number of vaccinated individuals supposed to not acquire sufficient immunity to become permanently recovered due to the short-term effect of vaccines. A basic reproduction number, when the control is defined as a constant parameter, is calculated. Since we also need to find the optimal values of this impulse control when it is defined as a function of time, we start by stating a general form of an impulse version of Pontryagin’s maximum principle that can be adapted to our case, and then we apply it to our model. Finally, we provide our numerical simulations that are obtained via an impulse progressive-regressive iterative scheme with fixed intervals between impulse times (theoretical example of an impulse at each week), and we conclude with a discussion of our results.

1. Introduction

In 1998, Shulgin et al. [1], highlighted the potential of pulse vaccination strategy in eradicating some epidemics at relatively low values of vaccination compared to conventional strategies as theoretically proven also in [2]. A theoretical examination of these results has been published after in [3]. In the first mentioned reference and more recently in paper of Yang and Xiao in [4], the authors discussed the effectiveness and importance of such vaccination regimes by referring to many theoretical and applicable results in the literature, namely the work of De Quadros et al. [5] in case of poliomyelitis, and work of Albert Bruce Sabin in [6] in cases of measles.

In 2001, Brauer et al. in [7], referred to the pulse vaccination suggested in [1,2,3] as a novel theoretical way to avoid recurring outbreaks of an epidemic. After the promising results presented in these references, many researchers showed more interest in the study of the complexity of epidemic models that involved this type of control measures, as in work of Alberto d’Onofrio in [8] where it was numerically observed through the use of SEIR model that pulse vaccination program is slightly more efficient than the traditional continuous one, Zhou and Liu in [9] with similar results in case of SIS model, Zeng and Chen in [10] where it is understood that SIRS model becomes much more complex under such vaccination policies while these authors concluded the same with Sun in [11] in case of SIR model. Extended studies with other different modeling frameworks have been published after, as in work of Gakkhar and Negi in [12] where the authors studied the complexity of SIRS system with non-monotonic incidence rate, Zhang and Teng in [13,14] where the two references results both indicated that a large pulse vaccination rate will help to eradicate the disease using SEIRS and SIRVS systems respectively, Li and Cui in [15] who compared between results of constant and pulse vaccination based on a comparison between pulse interval and a critical value of pulsing period, Zhao and Meng in [16] that treated the application of pulse vaccination to vertically transmitted diseases as considered earlier in other references as in [17,18], Pei et al. in [19] where it had been showed that the quarantine measure helps in eradicating the disease with similar conclusions as in [13,14] but through a delayed SEIQR model, Qiao et al. in [20] exhibiting the potential of pulse vaccination strategy in eliminating viral infectious diseases, see also other interesting literature in [21,22,23,24].

There have been intensive research papers about optimal vaccination strategies obtained through optimal control methods for many diseases prevention, see ref. [25] in case of dengue, ref. [26] for Influenza, refs. [27,28] for malaria, ref. [29] for Hepatitis B, ref. [30] for Ebola, ref. [31] for cholera. In these papers, the control is considered to be a continuous function in time, while the optimal control problem is resolved using classical theorems of sufficient and necessary conditions of optimality stated in [32,33] respectively. The optimization subject to such epidemic models, has mostly aimed to minimize the control policy cost associated with vaccination; however, as far as we know, there have been no studies that tried to resolve this problem in case of pulses with clearest and simplest methodology for better practicability. This means there is a need for a numerical method that applies exactly most findings of the theoretical calculus related to the necessary conditions of optimality with possibly clear analytical formulation of the sought optimal impulse control and without the help of any type of gradient methods. In 1978, Balquière provided in [34], an impulsive version for these optimality conditions, but there is still a lack in research for larger interpretation of these results, especially in the numerical part. Here, we try to develop an algorithm that will simplify more the resolution of these problems, often referred to as optimal impulsive control problems [35]. Further results about existence of solutions to these problems, have been discussed in [36,37].

We consider an optimal control problem subject to a state system that involves a finite number of impulse times in which an impulse control function is applied. More concretely, the control will represent the vaccination ratio as a variable in function of time and that is applied to a proportion susceptible to infection, while its optimal values will be sought at each time of vaccine application for better management of the pulse vaccination strategy.

To reach our purpose, we consider a susceptible-infected-removed (SIR) epidemic model with an additional compartment to represent the number of the vaccinated population. The proposed model is in the form of a differential system which contains a parameter that characterizes an impulse vaccination strategy followed for preventing the spread of an epidemic.

Since not every vaccine is necessarily $100\%$ effective as explained in [38], we choose to study dynamics of the considered population in this case. Then, an optimal impulse control strategy is suggested for this model using an impulse version of Pontryagin’s maximum principle. Finally, impulse progressive-regressive iterative schemes are used for finding the results.

2. Epidemic Control Model With Short-Term Immunity

2.1. Model Description in Continuous-Time Case and with Constant Control

In this part, we consider an epidemic control model where the vaccination strategy is described using a constant control parameter with the hypothesis that the used vaccines have a limited effect on people who receive it.

We need first to define the following four compartments

$\dot{S}$:

the number of people susceptible to infection or who are not yet infected,

$\dot{V}$:

the number of susceptible people who are temporary vaccinated so they cannot move to the removed class due to the limited effect of vaccine [38,39,40],

$\dot{I}$:

the number of infected people who are capable of spreading the epidemic to those in the susceptible and temporary controlled categories,

$\dot{R}$:

the number of removed people from the epidemic.

For all t belonging to an interval $[0,T]$ and when the vaccination policy is only presented by a constant parameter of control $\theta$, the model takes the form of the following differential system

$$\left\{\begin{array}{ccc}\dot{S}(t)\hfill & =\hfill & \mathsf{\Pi }(t)-\beta S(t)I(t)-a\theta S(t)-\mu S(t)\hfill \\ \dot{V}(t)\hfill & =\hfill & a\theta S(t)-b\beta V(t)I(t)-\mu V(t)\hfill \\ \dot{I}(t)\hfill & =\hfill & \beta (S(t)+bV(t))I(t)-\gamma I(t)-\mu I(t)\hfill \\ \dot{R}(t)\hfill & =\hfill & \gamma I(t)-\mu R(t)\hfill \end{array}\right.$$

(1)

with initial conditions $S(0)=S_0>0$, $V(0)=V_0\ge 0$, $I(0)=I_0\ge 0$ and $R(0)=R_0\ge 0$ and where $\mathsf{\Pi }(t)=\mu N(t)$ with $N(t)=S(t)+V(t)+I(t)+R(t)$, giving the newborn people.

$a\theta$$(0\le a\le 1)$ is the fraction of controlled people with “a” modeling the reduced chances of a susceptible individual to be controlled.

Also, $\beta =\frac{\kappa }{N(t)}$ with $\kappa$ the infection transmission rate, $\mu$ is the natural death rate, $b\beta$ $(0\le b\le 1)$ is the infection transmission rate even in the presence of $\theta$ with “b” modeling the reduced chances of a temporary controlled individual to be infected, and $\gamma$ is the recovery rate.

The population size N is constant because $\dot{N}(t)=\dot{S}(t)+\dot{V}(t)+\dot{I}(t)+\dot{R}(t)=0$, hence, $\mathsf{\Pi }(t)=\mathsf{\Pi }=\mu N_0$ knowing that $N(0)=N_0$.

We note that $\theta =0$ will denote no vaccination and $\theta =1$ will denote the use of vaccination with an initial number of susceptible people $S(0)$.

In [38], the authors devised a discrete-time version of (1) under the reason that in epidemics, data are collected at discrete times, while in [39,40], the model has been limited to its following continuous-time version

$$\left\{\begin{array}{ccc}\dot{S}(t)\hfill & =\hfill & \mathsf{\Pi }(t)-\beta S(t)I(t)-a\theta (t)S(t)-\mu S(t)\hfill \\ \dot{V}(t)\hfill & =\hfill & a\theta (t)S(t)-b\beta V(t)I(t)-\mu V(t)\hfill \\ \dot{I}(t)\hfill & =\hfill & \beta (S(t)+bV(t))I(t)-\gamma I(t)-\mu I(t)\hfill \\ \dot{R}(t)\hfill & =\hfill & \gamma I(t)-\mu R(t)\hfill \end{array}\right.$$

(2)

However, since vaccinations are often neither given continuously, nor at daily or very close discrete times, an impulse version of such models, seems more realistic. This leads us to think about devising an epidemic model that takes into account such considerations while providing theoretical and numerical steps for the resolution of an impulse optimal control problem.

Let first present the model (1) when there is use of pulse vaccinations.

2.2. Model Formulation in Impulsive Case with Control Function

Let $t_1$, $t_2$, …, $t_n$ be the times the vaccinations are applied. n is their total number.

Since our goal in this paper is to find the optimal values of pulse vaccination rate $\theta$ at $t_{i}i=1,\dots ,n$, then, we start with converting the continuous-time control $\theta (t)$ into the discrete summation ${\sum }_{i=1}^n\delta (t-t_i)\theta (t_i)$. The continuous-time model in (2) becomes defined by the following system with impulse vaccinations

$$\left\{\begin{array}{ccc}\dot{S}(t)\hfill & =\hfill & {\displaystyle \mathsf{\Pi }-\beta S(t)I(t)-a\sum _{i=1}^n\delta (t-t_i)\theta (t_i)S(t)-\mu S(t)}\hfill \\ \dot{V}(t)\hfill & =\hfill & {\displaystyle a\sum _{i=1}^n\delta (t-t_i)\theta (t)S(t)-b\beta V(t)I(t)-\mu V(t)}\hfill \\ \dot{I}(t)\hfill & =\hfill & \beta (S(t)+bV(t))I(t)-\gamma I(t)-\mu I(t)\hfill \\ \dot{R}(t)\hfill & =\hfill & \gamma I(t)-\mu R(t)\hfill \end{array}\right.$$

(3)

where $\delta (t-t_i)$ is a type of Dirac-delta function and each term of the summation above realizes the following property

$$\delta (t-t_i)\theta (t_i)=\left\{\begin{array}{cc}\theta (t_i)\hfill & \mathrm{if}t=t_i\hfill \\ 0\hfill & \mathrm{if}t\ne t_i\hfill \end{array}\right.$$

(4)

In fact, the summation ${\sum }_{i=1}^n\delta (t-t_i)\theta (t_i)$ in (3) can be considered to be a sum of weighted Kronecker-delta functions as the term (4) is a shifted delta function with a specified weight. More generally, the summation as defined here, refers to a bounded discrete signal and the delta function is defined in this case by

$$\delta (t-t_i)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}t=t_i\hfill \\ 0\hfill & \mathrm{if}t\ne t_i\hfill \end{array}\right.or\delta (t)=\left\{\begin{array}{cc}1\hfill & \mathrm{if}t=0\hfill \\ 0\hfill & \mathrm{if}t\ne 0\hfill \end{array}\right.$$

(5)

while an integral form refers to real continuous-time case and the delta function is defined as

$$\delta (t)=\left\{\begin{array}{cc}+\infty \hfill & \mathrm{if}t=0\hfill \\ 0\hfill & \mathrm{if}t\ne 0\hfill \end{array}\right.$$

(6)

with the famous properties ${\int }_{-\infty }^{+\infty }\delta (t)dt=1$, ${\int }_{-\infty }^{+\infty }\theta (t)\delta (t)dt=\theta (0)$ and ${\int }_{-\infty }^{+\infty }\theta (t-\alpha )\delta (t)=\theta (\alpha ),$ $\alpha \in \mathbb{R}$, but this is not the case here.

Finally, and contrary to the form of control in (2), $\theta (t)$ will be defined as the whole sequence in one single function ${\sum }_{i=1}^n\delta (t-t_i)\theta (t_i)$ and when seeking its optimal values, we will need the values of $\theta (t_i)$.

Thus, we understand the following conditions.

If $t\ne t_i$ $\forall i$, then the system (3) is changed to

$$\left\{\begin{array}{ccc}\dot{S}(t)\hfill & =\hfill & \mathsf{\Pi }-\beta S(t)I(t)-\mu S(t)\hfill \\ \dot{V}(t)\hfill & =\hfill & -b\beta V(t)I(t)-\mu V(t)\hfill \\ \dot{I}(t)\hfill & =\hfill & \beta (S(t)+bV(t))I(t)-\gamma I(t)-\mu I(t)\hfill \\ \dot{R}(t)\hfill & =\hfill & \gamma I(t)-\mu R(t)\hfill \end{array}\right.$$

(7)

and if $t=t_i$ $\forall i$, then we can obtain the following equations

$$\left\{\begin{array}{ccc}S(t_i^{+})\hfill & =\hfill & S(t_i)-a\theta (t_i)S(t_i)\hfill \\ V(t_i^{+})\hfill & =\hfill & V(t_i)+a\theta (t_i)S(t_i)\hfill \\ I(t_i^{+})\hfill & =\hfill & I(t_i)\hfill \\ R(t_i^{+})\hfill & =\hfill & R(t_i)\hfill \end{array}\right.$$

(8)

where $t_i$ representing the vaccination time and $\theta (t_i)$ the impulse control at this time, while $S(t_i)=S(t_i^{-})$, $V(t_i)=V(t_i^{-})$, $I(t_i)=I(t_i^{-})$ and $R(t_i)=R(t_i^{-})$ are defined as the old numbers of susceptible, vaccinated, infected, and removed people respectively, while $S(t_i^{+})$, $V(t_i^{+})$, $I(t_i^{+})$ and $R(t_i^{+})$ defining these same numbers just after $t_i$. Thus, (8) means we have abandoned the old vaccination and used a new vaccination for the number of susceptible people $S(t_i)$.

Before the presentation of our optimal control approach, let provide first a result about the stability of these systems when the control parameter is still defined as a constant rate.

3. Basic Reproduction Number

In system (7) and (8), since the three first equations are independent of the last equation, the study of stability can be reduced and done for the following system

$$\left\{\begin{array}{ccc}\dot{S}(t)\hfill & =\hfill & \mathsf{\Pi }-\beta S(t)I(t)-\mu S(t)\hfill \\ \dot{V}(t)\hfill & =\hfill & -b\beta V(t)I(t)-\mu V(t)\hfill \\ \dot{I}(t)\hfill & =\hfill & \beta (S(t)+bV(t))I(t)-\gamma I(t)-\mu I(t)\hfill \\ S(t_i^{+})\hfill & =\hfill & S(t_i)-a\theta S(t_i)\hfill \\ V(t_i^{+})\hfill & =\hfill & V(t_i)+a\theta S(t_i)\hfill \\ I(t_i^{+})\hfill & =\hfill & I(t_i)\hfill \end{array}\right.$$

(9)

which has a unique and piecewise continuous solution in $(t_i,t_{i+1}]$ [41,42].

If there is no infection, i.e., $I(0^{+})$ then $I(t)\equiv 0,t\ge 0$, and system (9) becomes

$$\left\{\begin{array}{ccc}\dot{S}(t)\hfill & =\hfill & \mathsf{\Pi }-\mu S(t)\hfill \\ \dot{V}(t)\hfill & =\hfill & -\mu V(t)\hfill \\ S(t_i^{+})\hfill & =\hfill & S(t_i)-a\theta S(t_i)\hfill \\ V(t_i^{+})\hfill & =\hfill & V(t_i)+a\theta S(t_i)\hfill \end{array}\right.$$

(10)

whose solution on $(0,h]$ with h, the time between two consecutive impulse vaccinations, is

$$\left\{\begin{array}{ccc}S(t)\hfill & =\hfill & N_0+({\tilde{S}}_0-N_0)e^{-\mu t}\hfill \\ V(t)\hfill & =\hfill & {\tilde{V}}_0{e}^{-\mu t}\hfill \end{array}\right.$$

(11)

with ${\tilde{S}}_0=S(0^{+})$ and ${\tilde{V}}_0=V(0^{+})$.

In fact, the first equation in (10) implies that $\dot{S}(t)=\mu N_0-\mu S(t)=\mu (N_0-S(t))$ as noted earlier in the definition of $\mathsf{\Pi }$ in Section 2. By passage to the integral with letting just $\mu$ in the right-hand side, we obtain, $-ln(N_0-S(t))=\mu t+c$, then, $N_0-S(t)=c{e}^{-\mu t}$, which gives $S(t)=N_0-c{e}^{-\mu t}$. At ${\tilde{S}}_0$, we have ${\tilde{S}}_0=N_0-c$, then, $c=N_0-{\tilde{S}}_0$ and by replacing it in $S(t)$, we obtain the result. As for the second equation, this is simpler after integration and we can find that constant c is equal to ${\tilde{V}}_0$.

Let $S(h^{+})=S(0^{+})$ and $V(h^{+})=V(0^{+})$, we can solve ${\tilde{S}}_0$, ${\tilde{V}}_0$ and find a unique positive periodic solution $(\tilde{S}(t),\tilde{V}(t))$ of system (10).

$$\left\{\begin{array}{ccc}\tilde{S}(t)\hfill & =\hfill & N_0+({\tilde{S}}_0-N_0)e^{-\mu (t-t_i)},t\in (t_i,t_{i+1}]\hfill \\ \tilde{V}(t)\hfill & =\hfill & {\tilde{V}}_0{e}^{-\mu (t-t_i)},t\in (t_i,t_{i+1}]\hfill \end{array}\right.$$

(12)

where

$$\left\{\begin{array}{ccc}{\tilde{S}}_0\hfill & =\hfill & {\displaystyle \frac{N_0(1-a\theta )(1-e^{-\mu h})}{1-(1-a\theta )e^{-\mu h}}}\hfill \\ {\tilde{V}}_0\hfill & =\hfill & {\displaystyle \frac{a\theta (N_0+({\tilde{S}}_0-N_0)e^{-\mu h}}{(1-e^{-\mu h})}}\hfill \end{array}\right.$$

(13)

Then, we find the basic reproduction number associated with system (9) as

$${\displaystyle R_0^I=\frac{1}{h}{\int }_0^h\frac{\beta (\tilde{S}(t)+b\tilde{V}(t))}{\mu +\gamma }dt}$$

(14)

Based on similar results in [43,44,45,46,47], we announce the following theorem

Theorem 1.

The disease free periodic solution $(\tilde{S}(t),\tilde{V}(t),0)$ is globally asymptotically stable if $R_0^I<1$ and the disease is permanent if $R_0^I>1$ meaning the existence of a positive constant $\delta >0$ such that any solution $(S(t),V(t),I(t))$ of system (9) with initial values $S(0^{+})\ge 0,V(0^{+})\ge 0,I(0^{+})\ge 0$ satisfies ${lim}_{t\to \infty }infI(t)\ge \delta$.

This theorem shows that system (9) exhibits similar threshold dynamics to the continuous-time case treated in [39].

4. Optimal Impulse Control Approach

4.1. Theoretical Framework

In this part of paper, we define our objective by determining the optimal values of the impulse control $\theta (t_i)$ that minimizes a given criterion, but before this, we present hereafter the theoretical steps for reaching this goal.

With a combination of results in [48,49], we try here to derive an analogous necessary conditions of optimality that will be associated with our particular form of the impulsive system (3) which does not contain any continuous-time control as usually taken into consideration in general theoretical formulation in the literature, but only an impulse control here.

We start then first by presenting a general formulation of our optimal control problem that can be adapted to the form (7) and (8) as follows

$${\displaystyle \underset{n,t_i,\theta (t_i)}{max}{\int }_0^{T}F(t,x(t))+\sum _{i=1}^{n}G(x(t_i^{-},\theta (t_i),t_i))+\varphi (x(T^{+}))}$$

(15)

subject to

$$\dot{x}(t)=f(t,x(t)),t\notin \{t_1,\dots ,t_n\}$$

$$x(t_i^{+})-x(t_i^{-})=g(x(t_i^{-},\theta (t_i),t_i)),i\in \{1,\dots ,n\}$$

$$x\in {\mathbb{R}}^N,\theta (t_i)\in \mathsf{\Theta },x(0^{-})=x_0,t_i\in [0,T]$$

where x is the state variable, a piecewise continuous function of time, and $\theta (t_i)$ is the impulse control variable. n is the number of pulses, $t_i$ is the instant of the ith pulse and $t_i^{-}$, $t_i^{+}$ are the instants just before and after the impulse, (i.e., $x(t_i^{-})$, $x(t_i^{+})$ represent the first and second limit sides of x respectively). The final time is noted $T>0$ while $T^{+}$ is the instant that comes just after T.

The system gain is given by $F(t,x(t))$, $G(x(t_i^{-},\theta (t_i),t_i))$ is the gain function associated with the ith pulse, and $\varphi (x(T^{+}))$ is the final cost function associated with the system just after T.

Finally, $f(t,x(t))$ is the continuous change of the state variable through time between pulses points and $g(x(t_i^{-},\theta (t_i),t_i))$ is the function that represents the instantaneous or finite change of the state variable when there is a pulse.

We admit that $\mathsf{\Theta }$ is a bounded convex control set, and we impose that $F,f,g$ and G are continuously differentiable functions in x on ${\mathbb{R}}^N$ and in $\theta (t_i)$ on $\mathsf{\Theta }$, $\varphi (x(T^{+}))$ is continuously differentiable in $x(T^{+})$ on ${\mathbb{R}}^n$,and that g and G are continuous on t. Finally, when there is no pulse, i.e., $\theta (t_i)=0$, we assume that $g(x(t),0,t)=0$ for all x and t.

As we have defined the general formulation of our objective that can be associated with the form of impulsive system (3), we follow similar steps as in [50,51,52] to obtain the following theorem which states the necessary conditions of optimality associated with our special case of optimal impulse control problem (15).

Theorem 2.

Let $(x^{*}(t),n,t_1^{*},\dots ,t_k^{*},{\theta }^{*}(t_1),\dots ,{\theta }^{*}(t_k))$ be optimal solutions of the maximization problem (15), then there exists an adjoint variable such that the following conditions are satisfied

$$\begin{array}{ccc}\hfill \dot{\lambda }(t)& =& {\displaystyle -\frac{\partial H}{\partial x}(t,x^{*}(t),\lambda (t)).}\hfill \end{array}$$

(16)

In the points of pulses, we have

$$H^I(t_i,x^{*}(t),{\theta }^{*}(t_i))\ge H^I(t_i,x(t_i),\theta )\forall \theta \in \mathsf{\Theta }$$

(17)

$$\begin{array}{ccc}\hfill \lambda (t_i^{*+})-\lambda (t_i^{*-})& =& {\displaystyle -\frac{\partial H^I}{\partial x}(t_i^{*},x^{*}(t_i^{*-}),{\theta }^{*}(t_i),\lambda (t_i^{*+}))}\hfill \end{array}$$

(18)

Proof. 

In [34], the author developed a maximum principle theorem for necessary conditions of optimality to resolve the problem (15) but with a continuous-time control function in the state equation when $t\ne t_i$. Here, since we have no continuous-time control in system (7), we proceed in our generalized form as follows.

First, we consider the Hamiltonian function defined by

$$H(t,x(t),\lambda (t))=F(t,x(t))+\lambda {(t)}^{T}f(t,x(t))$$

and the impulse Hamiltonian function defined by

$$H^I(t,x(t),\theta (t_i))=G(t,x(t),\theta (t_i))+\lambda (t^{+})g(t,x(t),\theta (t_i))$$

where $\lambda (t)$ represents the adjoint state variable.

Please note that (16) is equivalent to the existence of the adjoint system in the classical continuous-time version of maximum principle [33]. As for (17) and (18), it is directly obtained from results in [34] and the only difference here in (17) is that the sign ≤ is replaced by ≥ because (15) is presented as a maximization problem. An analogous condition for the Hamiltonian H is necessary if the state equation also includes a continuous-time control, otherwise, the only sought optimal control here is the impulse one that intervenes when $t=t_i$. In this remark, a function $H(t,x^{*}(t),\lambda (t))$ still has a meaning because instead of defining an optimal trajectory $x^{*}$ for both continuous and impulse controls, $x^{*}$ will be just associated now to ${\theta }^{*}$ only. This also implies that another condition may involve both classical and impulse Hamiltonian functions H and $H^I$ and which can be defined by

$$\begin{array}{c}\hfill {\displaystyle H(t_i,x^{*}(t_i^{+}),\lambda (t_i^{+}))+\frac{\partial H^I}{\partial t}(t_i,x^{*}(t_i^{+}),{\theta }^{*}(t_i^{+}))}\\ \hfill {\displaystyle =H(t_i,x^{*}(t_i),\lambda (t_i))+\frac{\partial H^I}{\partial t}(t_i,x^{*}(t_i),{\theta }^{*}(t_i))}\end{array}$$

(19)

Since we will deal with an autonomous system, the condition (19) is reduced to

$$\begin{array}{c}\hfill H(t_i,x^{*}(t_i^{+}),\lambda (t_i^{+}))=H(t_i,x^{*}(t_i),\lambda (t_i))\end{array}$$

and it means that H is continuous at all times where ${\theta }^{*}$ is applied [49]. □

4.2. Optimal Impulse Control Problem

The most evident objective in epidemic control as defined for models in [38,39], concerns the minimization of the number of infected people while minimizing the cost of vaccination. Thus, for our case, we choose to define our objective function to be minimized as

$${\displaystyle J(\theta )={\int }_0^{T}AI(t)dt+\sum _{i=1}^n\frac{B}{2}{\theta }^2(t_i)}$$

(20)

In other words, we are seeking an optimal control ${\theta }^{*}$ such that

$$J({\theta }^{*})=min\{J(\theta )/\theta \in \mathsf{\Theta }\}$$

(21)

subject to

$$\left\{\begin{array}{c}{\displaystyle \dot{S}(t)=\mathsf{\Pi }-\beta S(t)I(t)-a\sum _{i=1}^n\delta (t-t_i)\theta (t_i)S(t)-\mu S(t)}\hfill \\ {\displaystyle \dot{V}(t)=a\sum _{i=1}^n\delta (t-t_i)\theta (t)S(t)-b\beta V(t)I(t)-\mu V(t)}\hfill \\ \dot{I}(t)=\beta (S(t)+bV(t))I(t)-\gamma I(t)-\mu I(t)\hfill \\ \dot{R}(t)=\gamma I(t)-\mu R(t)\hfill \\ S_0,V_0,I_0,R_{0}given\hfill \\ 0\le \theta (t)\le 1\hfill \end{array}\right.$$

(22)

Or, our optimal pulse control problem can be stated as follows

$$\left\{\begin{array}{ll}& {\displaystyle \underset{\theta \in \mathsf{\Theta }}{min}\left\{{\displaystyle {\int }_0^{T}AI(t)dt+\sum _{i=1}^n\frac{B}{2}{\theta }^2(t_i)}\right\}}\hfill \\ & subjectto\hfill \\ & \dot{S}(t)=\mathsf{\Pi }-\beta S(t)I(t)-\mu S(t)\hfill \\ & \dot{V}(t)=-b\beta V(t)I(t)-\mu V(t)\hfill \\ & \dot{I}(t)=\beta (S(t)+bV(t))I(t)-\gamma I(t)-\mu I(t)\hfill \\ & \dot{R}(t)=\gamma I(t)-\mu R(t)\hfill \\ whent\ne t_i\hfill & \\ and\hfill & \\ & S(t_i^{+})=S(t_i)-a\theta (t_i)S(t_i)\hfill \\ & V(t_i^{+})=V(t_i)+a\theta (t_i)S(t_i)\hfill \\ & I(t_i^{+})=I(t_i)\hfill \\ & R(t_i^{+})=R(t_i)\hfill \\ whent=t_i\hfill & \\ & S_0,V_0,I_0,R_{0}given\hfill \\ & 0\le \theta (t)\le 1\hfill \end{array}\right.$$

To resolve problem (22), we follow these steps.

First, we construct the Hamiltonian H without control as the function defined as

$$\begin{array}{ll}& H(S(t),V(t),I(t),R(t),{\lambda }_1(t),{\lambda }_2(t),{\lambda }_3(t))\hfill \\ =& AI(t)+{\lambda }_1(t)(\mathsf{\Pi }-\beta S(t)I(t)-\mu S(t))-{\lambda }_2(t)(b\beta V(t)I(t)+\mu V(t))\hfill \\ & +{\lambda }_3(t)(\beta (S(t)+bV(t))I(t)-\gamma I(t)-\mu I(t))+{\lambda }_4(t)(\gamma I(t)-\mu R(t))\hfill \end{array}$$

(23)

and second, we construct the impulse Hamiltonian function $H^p$ defined as

$$\begin{array}{c}\hfill H^p(S(t),V(t),I(t),R(t),\theta (t))=\frac{B}{2}{\theta }^2(t)-a{\lambda }_1(t^{+})\theta (t)S(t)+a{\lambda }_2(t^{+})\theta (t)S(t)\end{array}$$

(24)

We prefer the notation $H^p$ to $H^I$ just to avoid any ambiguity especially when we will need the derivative of this function with respect to function of infection I.

We can announce the impulse maximum principle as follows.

Theorem 3.

Given an impulse optimal control ${\theta }^{*}(t_i),i=1,\dots ,n$ where $t_i>t_{i-1}>0$ and which minimizes (20) along with the optimal trajectories $S^{*}$, $V^{*}$, $I^{*}$ and $R^{*}$ associated with the differential system in (22), then there exist adjoint variables ${\lambda }_k$, $k=1,2,3,4$ as notations of ${\lambda }_k(t)$ and which satisfy for $t\ne t_i$ the following adjoint differential system

$$\left\{\begin{array}{ccc}\dot{{\lambda }_1}(t)\hfill & =\hfill & {\lambda }_1(t)(\beta I^{*}(t)+\mu )-\beta {\lambda }_3(t)I^{*}(t)\hfill \\ \dot{{\lambda }_2}(t)\hfill & =\hfill & {\lambda }_2(t)(b\beta I^{*}(t)+\mu )-b{\lambda }_3(t)\beta I^{*}(t)\hfill \\ \dot{{\lambda }_3}(t)\hfill & =\hfill & -A+{\lambda }_1(t)\beta S^{*}(t)+b{\lambda }_2(t)\beta V^{*}(t)-{\lambda }_3(t)(\beta (S^{*}(t)+b{V}^{*}(t))-\mu -\gamma )-{\lambda }_4(t)\gamma \hfill \\ \dot{{\lambda }_4}(t)\hfill & =\hfill & {\lambda }_4(t)\mu \hfill \end{array}\right.$$

(25)

with the transversality conditions ${\lambda }_k(T)=0$, $k=1,2,3,4$ and we have,

$$\left\{\begin{array}{ccc}{\lambda }_1(t_i)\hfill & =\hfill & {\lambda }_1(t_i^{+})-a{\lambda }_1(t_i^{+})\theta (t_i)+a{\lambda }_2(t_i^{+})\theta (t_i)\hfill \\ {\lambda }_2(t_i)\hfill & =\hfill & {\lambda }_2(t_i^{+})\hfill \\ {\lambda }_3(t_i)\hfill & =\hfill & {\lambda }_3(t_i^{+})\hfill \\ {\lambda }_4(t_i)\hfill & =\hfill & {\lambda }_4(t_i^{+})\hfill \end{array}\right.$$

(26)

Furthermore, we have

$${\theta }^{*}(t_i)=min\left(max\left({\displaystyle 0,\frac{a{S}^{*}(t_i)({\lambda }_1(t_i^{+})-{\lambda }_2(t_i^{+}))}{B}}\right),1\right)$$

(27)

Proof. 

We have the Hamiltonian function H is defined in (23) as

$$\begin{array}{c}H(S(t),V(t),I(t),R(t),{\lambda }_1(t),{\lambda }_2(t),{\lambda }_3(t))\hfill \\ =AI(t)+{\lambda }_1(t)(\mathsf{\Pi }-\beta S(t)I(t)-\lambda S(t))-{\lambda }_2(t)(b\beta V(t)I(t)+\lambda V(t))\hfill \\ +{\lambda }_3(t)(\beta (S(t)+bV(t))I(t)-\gamma I(t)-\lambda I(t))+{\lambda }_4(t)(\gamma I(t)-\lambda R(t))\hfill \end{array}$$

(28)

Then, using the impulse version of Pontryagin’s maximum principle in [48], then we have, for $t\ne t_i$

$$\begin{array}{ccc}\hfill \dot{{\lambda }_1}(t)& =& {\displaystyle -\frac{\partial H}{\partial S}(S(t),V(t),I(t),R(t),{\lambda }_1(t),{\lambda }_2(t),{\lambda }_3(t))}\hfill \\ & =& {\lambda }_1(t)(\beta I^{*}(t)+\mu )-\beta {\lambda }_3(t)I^{*}(t)\hfill \\ \hfill \dot{{\lambda }_2}(t)& =& {\displaystyle -\frac{\partial H}{\partial V}(S(t),V(t),I(t),R(t),{\lambda }_1(t),{\lambda }_2(t),{\lambda }_3(t))}\hfill \\ & =& {\lambda }_2(t)(b\beta I^{*}(t)+\mu )-b{\lambda }_3(t)\beta I^{*}(t)\hfill \\ \hfill \dot{{\lambda }_3}(t)& =& {\displaystyle -\frac{\partial H}{\partial I}(S(t),V(t),I(t),R(t),{\lambda }_1(t),{\lambda }_2(t),{\lambda }_3(t))}\hfill \\ & =& -A+{\lambda }_1(t)\beta S^{*}(t)+b{\lambda }_2(t)\beta V^{*}(t)-{\lambda }_3(t)(\beta (S^{*}(t)+b{V}^{*}(t))-\mu -\gamma )-{\lambda }_4(t)\gamma \hfill \\ \hfill \dot{{\lambda }_4}(t)& =& {\displaystyle -\frac{\partial H}{\partial R}(S(t),V(t),I(t),R(t),{\lambda }_1(t),{\lambda }_2(t),{\lambda }_3(t))}\hfill \\ & =& {\lambda }_4(t)\mu \hfill \end{array}$$

while the transversality conditions defined as minus the derivative of the final gain function with respect to the state variables S, $C_S$, I, and R. Since the final gain function in (20) does not contain any term of these variables, then ${\lambda }_k(T)=0$, $k=1,2,3,4$.

Since the impulse Hamiltonian function $H^p$ is defined in (24) as

$$\begin{array}{c}\hfill H^p(S,V,I,R,\theta )=\frac{B}{2}{\theta }^2-a{\lambda }_1(t^{+})\theta S+a{\lambda }_2(t^{+})\theta S\end{array}$$

Then, we have also at $t=t_i$

$$\left\{\begin{array}{ccc}{\lambda }_1(t_i)\hfill & =\hfill & {\lambda }_1(t_i^{+})+H_S^p(S,V,I,R,\theta )={\lambda }_1(t_i^{+})-a{\lambda }_1(t_i^{+})\theta (t_i)+a{\lambda }_2(t_i^{+})\theta (t_i)\hfill \\ {\lambda }_2(t_i)\hfill & =\hfill & {\lambda }_2(t_i^{+})+H_V^p(S,V,I,R,\theta )={\lambda }_2(t_i^{+})\hfill \\ {\lambda }_3(t_i)\hfill & =\hfill & {\lambda }_3(t_i^{+})+H_I^p(S,V,I,R,\theta )={\lambda }_3(t_i^{+})\hfill \\ {\lambda }_4(t_i)\hfill & =\hfill & {\lambda }_4(t_i^{+})+H_R^p(S,V,I,R,\theta )={\lambda }_4(t_i^{+})\hfill \end{array}\right.$$

The optimality condition at $\theta ={\theta }^{*}$ implies that ${\displaystyle \frac{\partial H^p}{\partial \theta (t_i)}=0}$.

Then, after setting $S=S^{*}$, we have

$${\displaystyle B\theta (t_i)-aS(t_i){\lambda }_1(t_i^{+})+aS(t_i){\lambda }_2(t_i^{+})=0\Rightarrow \theta (t_i)=\frac{aS(t_i)({\lambda }_1(t_i^{+})-{\lambda }_2(t_i^{+}))}{B}.}$$

By using the bounds of control defined in (22) and similar detailed calculus as in [38,39], we obtain the following characterization

$${\theta }^{*}(t_i)=min\left(max\left({\displaystyle 0,\frac{a{S}^{*}(t_i)({\lambda }_1(t_i^{+})-{\lambda }_2(t_i^{+}))}{B}}\right),1\right)$$

 □

5. Numerical Results

As mentioned in introduction, we aim in this numerical part, to propose a numerical method that applies exactly the findings of the last theorem. Since all necessary conditions of optimality along with the characterization of the optimal control, are all found explicitly in this theorem, we do not need any type of gradient method.

In fact, the four steps of numerical calculus associated with the resolution of our impulse optimal control problem (22)–(25), are described in Algorithm 1.

Figure 1 and Figure 2 depict $SVIR$ dynamics in the absence and presence of control. As for Figure 3, we exhibit the recommended pulses.

When we talk about the case when there is no vaccination strategy, we mean that we deal with the following continuous-time SIR differential equation

$$\left\{\begin{array}{ccc}\dot{S}(t)\hfill & =\hfill & \mathsf{\Pi }-\beta S(t)I(t)-\mu S(t)\hfill \\ \dot{I}(t)\hfill & =\hfill & \beta S(t)I(t)-\gamma I(t)-\mu I(t)\hfill \\ \dot{R}(t)\hfill & =\hfill & \gamma I(t)-\mu R(t)\hfill \end{array}\right.$$

(29)

for all time $t\in [0,T]$, and this system clearly describes the continuous spread of the epidemic in a population without any control intervention. The behaviors of $S(t),I(t)$ and $R(t)$ for all time $t\in [0,T]$ in this case, are presented with the blue curves in the three figures below. As for the case when we follow an impulse vaccination strategy subject to (22), the behaviors of $S^{*}(t),I^{*}(t)$ and $R^{*}(t)$ for all time $t\in [0,T]$, are presented with the red curves.

Algorithm 1: Resolution steps of the two-point boundary value optimal control problem (22)–(25)

Step 0:  Guess an initial estimation to $\theta$. Step 1:  Use the initial condition $S(0)$, $V(0)$, $I(0)$ and $R(0)$ and the stocked values by $\theta$.  Find the optimal states $S^{*}$, $V^{*}$, $I^{*}$ and $R^{*}$ which iterate forward in (23) in the two cases when $t\ne t_i$ and when $t=t_i$. A progressive fourth order Runge-Kutta scheme is incorporated for finding the state solutions from 0 to a fixed number of iterations m. Step 2:  Use the stocked values by $\theta$ and the transversality conditions (26).  Find the adjoint variables ${\lambda }_k$ for $k=1,2,3,4$ which iterate backward in (25) when $t\ne t_i$ and in (26) when $t=t_i$. A regressive fourth order Runge-Kutta scheme is incorporated for finding the costate solutions from m to 0. Step 3:  Update the control using new S, V, I, R and ${\lambda }_k$ for $k=1,2,3,4$ in the characterization of ${\theta }^{*}$ as presented in(27) Step 4:  Test the convergence based on the forward-backward sweep method criterion. If the values of the sought variables in this iteration and the final iteration are sufficiently small, check out the recent values as solutions. If the values are not small, go back to Step 1.

From most epidemiological data we have revised until now, especially in the case of controllable epidemics that appeared at specific seasons, we concluded that it is preferable to avoid parameters values that can lead to claiming the occurrence of an exponential increase of infection and deny the exponential decrease that just comes after the infection reaches its maximal peak, we also prefer to not analyze the advantage of our approach on fractions, but instead on numbers of individuals. The parameters values that will be considered here, may certainly help to understand this phenomenon. In fact, there was always an urgent control measure which aimed to prevent the spread of epidemics with the best available resources. The problem as mentioned in [7], is that we cannot be sure that the infection will not appear again since in practice, it is nearly impossible to vaccinate most susceptible individuals of a population because not all people would come to be vaccinated and not all vaccinations are successful [7,38,39,40]. In the same context, El Kihal et al. discussed other possible reasons such as the constraint of limited resources as in [38] or infodemics as in [53]. Thus, a better way to avoid the recurrence of a disease outbreaks, is to follow a control strategy based on pulsed vaccinations.

From Figure 1a, we can observe that the susceptible number S decreases exponentially towards 37 individuals and this is mostly due to infection plus natural death as it can be understood from the first equation in (29). In fact, the number of infectives I in Figure 2a increases exponentially towards an important number of more than 54 individuals in during the first week and decreases thereafter to just 25 individuals in the sixth week because of the natural death and natural recovery as it can also be deduced from second equation of (29), but this also may be due to another type of control intervention missed here but not sufficient as infection remains important and needs to be controlled.

In Figure 1a, the optimal state $S^{*}$ decreases rapidly towards 38 individuals in the first week and this value is smaller than the one taken by S because of the optimal impulse control intervention. However, once we go more forward in time, the effect of control becomes negligible and this causes a slight increase of $S_{*}$ between 10 and 12 days, and then an simultaneous increase of the function $I^{*}$ in Figure 2a in the same period to 57 individuals in the first week but decreases again to 27 people and this implies to some small increases of $S^{*}$ compared to S, and we remark once there is a vaccination interruption just after the time of impulse, $I^{*}$ does a slight jump again and again causing repeated small decreases of $S^{*}$. Simultaneously, the number of removed people R in Figure 1b increases to 57 people before the third week because of the natural recovery or any effect of other missed control policy in removing the infectives, while $R^{*}$ reaches a value that is higher than this number with a maximal peak equaling to 62 in the third week.

From Figure 3, we can show that vaccination is recommended in all weeks except on the third time of impulsion. The optimal values of controls take the value 1 as maximal peak in these weeks while it is zero in week 3. In fact, as we can observe from Figure 1b and Figure 2a, since R and I functions do not oscillate during days that are $<14$ and $>28$ but they increase and decrease in this period respectively, then, there is no need to introduce any control in the third week.

The main result we can recapitulate from these first three figures is that all $S^{*},R^{*}$ and $I^{*}$ oscillated in all times of pulsed vaccinations, and once a vaccination is stopped, the infection goes up again in repeatedly increases with differentiated peaks, and respectively, the level of recovery goes down many times. As for $V^{*}$ in Figure 2b, it is very small because of the small values $S^{*}$ has taken in most times of vaccination.

6. Discussion and Conclusions

The introduction of a continuous-time control in system (29) as done in [39], would seem less important than the use of an impulse control when it comes to subject of vaccination since we know that vaccines are injected on a discrete basis. In fact, with an optimal impulse control strategy, we may succeed to obtain a reduction in the number of infected people (resp. an augmentation in the number of removed people) without following a control policy at all continuous-time t. From Figure 4, we remark there are no oscillations in the curves of I and R functions since there is no interruption in the control along $[0,T]$. Although less realistic, the continuous-time control also shows a positive effect on shape of R (resp. on shape of I) as the recovery has stabilized (resp. the infection has decreased) to small values a little before than the cases in Figure 1b and Figure 2a. However, vaccination managers would surely choose recommendations of Figure 4 since the main goal of the minimization problem (21) with respect to the definition in (20) has been reached in similar way with even many vaccination rest periods, and this also meets the conclusions of [8,9] that impulses are slightly better than administrating vaccines continuously. This is to say there is no need to have an optimal control in its maximal bound suggested overall t as observed in Figure 4. Talking about the bound 1 as maximal suggested value for the control either in Figure 4 or Figure 3, we mention this is due to the number of infected people that despite the decrease, its minimal value is still around 40. Thus, we try to resolve and discuss this problem in the two following figures.

In Figure 5 and Figure 6, we discuss the behaviors of $I(t)$ function and the recommended impulse controls when the infection transmission rate β varies. In Figure 5, when $\beta =0.002$, we can observe that the number of infected people decreases now to 20 individuals as minimal value, and simultaneously, the value of the optimal impulse control in the last week is reduced to 0.75. As for the case when $\beta =0.001$, $I(t)$ decreases to only 5 individuals while the value of the optimal control takes only 0.14 in the last week. In the last case of the same figure, we could now say that infection could be eradicated as long as we move forward in time since $I(t)$ is equal to approximately zero in the 49th day even no impulse control is suggested in this time but only a value of 0.55 in the previous week is sufficient, without forgetting the ones values of other controls in the first weeks.

As we wonder whether there is a possibility to avoid more impulses for lesser costs, we found in Figure 6 that by affecting smaller values to β than in the previous figure, this may be achieved. In fact, when $\beta =0.00006$, the infection reaches the zero value in the 35th day while there is no suggested control in the last and sixth weeks with negligible value of the control in week 5 and very small value in week 4, and with the value 1 as maximal peak in the first week only and the value 0.64 in the second one. Similar results are obtained when $\beta =0.00004$ with some reduction of 0.1 and 0.2 in the two first weeks, respectively. As for the last case, when $\beta =0.00001$, the two first impulse optimal controls take only 0.22 and 0.08 with no recommendations for other weeks. Thus, we deduce that the more the infection transmission rate is important, the fewer impulse controls are required to introduce and the infection reaches the zero value in minimal time.

In Figure 7, we discuss the behaviors of $I(t)$ function and the suggested impulse controls when the recovery rate γ varies. When $\gamma =0.2$, the infection decreases towards the zero value after only 20 days with the requirement of only two impulse controls taking the values 1 and 0.5 as maximal peaks in the two first weeks, respectively. Now, when $\gamma =0.3$, the infection lasts only 14 days and disappears with the introduction of only two impulse controls taking the values 0.48 and 0.07 in weeks 1 and 2, respectively. In the last case, when $\gamma =0.4$, the infection is eradicated after only 11 days with the effect of only one control in the first week and that takes a very small rate of 0.15 while a second control appears negligible in week 2. Thus, we deduce that the more the recovery rate is important, the fewer impulse controls are recommended while $I(t)$ function may reach 0 even before the second week of intervention.

From all the results exhibited in the last three figures with different values of γ and β, we can now agree with the conclusions of [13,14,19] that large pulse vaccination rate will help to eradicate the disease, but this remains correct only if the transmission rate is not very important since in the example of Figure 1, despite the distinguishable decrease of infection with respect to the case without control, the number of infection remains important even if the optimal control in Figure 4 takes its maximal bound.

Finally, the optimal control strategy subject to our SVIR hybrid epidemic model, showed promising results when applied for first time in this form to the problem of pulsed vaccinations. Such optimization approaches can be extended to other applications with better description of policies that are not followed on continuous basis but rather introduced just between fixed time intervals. In future work, the potential of this impulse control strategy may also be discussed in case we focus more on most promiscuous individuals only, as deduced through targeted immunization approach suggested in [54]. Another interesting work that can meet the idea of impulse vaccination, is in [55] where the authors found that the influential nodes are not vaccinated on continuous basis but just repeatedly.