Optimal Vaccination Campaigns Using Stochastic SIR Model and Multiobjective Impulsive Control

ABSTRACT

A multiobjective impulsive control scheme is proposed to answer how optimal vaccination campaigns should be implemented, regarding two conflicting targets: making the total number of infecteds small and the vaccination campaign as handy as possible. In this paper, a stochastic SIR model is used to better depict the characteristics of a disease in practical terms, where little influences may lead to sudden and unpredictable changes in the behavior of transmissions. This model is extended to analyze the effects of impulsive vaccinations in two phases: the transient regime control, taking into account the necessity to reduce the number of infected individuals to an acceptable level in a finite time, and the permanent regime control, that will act with fixed vaccinations to avoid another outbreak. A parallel version of NSGA-II is used as the multiobjective optimization machinery, considering both the probability of eradication and the vaccination campaign costs. The final result using the proposed framework shows nondominated policies that can guide public managers to decide which is the best procedure to be taken depending on the present situation.

Keywords:
planning of vaccination campaigns; multiobjective optimization; impulsive control; stochastic SIR

1 INTRODUCTION

One of the most studied models describing the behavior of a disease in a population is the endemic SIR (Susceptible-Infectious-Recovered) ¹⁹. It is based on deterministic nonlinear differential equations describing the influence of changes in one group on the other ones, in the same species, wherein the whole population is partitioned into three groups of Susceptible (individuals who are not infected, but may become), Infected (individuals who are sick and can transmit the disease to susceptible individuals) and Recovered (individuals who had been sick or vaccinated and are now immune to the contagion).

In previous works, some steps have been sketched toward the multiobjective optimization design of impulsive vaccination control of epidemics in a given time horizon. The work ⁷ proposed using non-fixed pulse vaccination in a deterministic SIR, via open-loop dynamic optimization methods. In a decision process, this can help to choose the specific policy to be implemented. Reference ¹¹ proposed a multiobjective optimization approach, minimizing both the infected individuals and the cost with vaccination. The nondominated solutions are validated in an Individual Based Model (IBM) to study a set of interval of confidence for each optimal policy strategy. Finally, the article ¹⁰ proposed an optimization scheme, implemented with a local search procedure, and a hash table is then adopted to cope with the optimization machinery in order to prevent variability. It is also proposed a control in two parts: choosing a set of optimal permanent regimes coupled with a set of temporary regimes. The validation of campaigns is made using a stochastic discrete-time model.

Evolutionary algorithms have become an important computational tool for dealing with objective functions which can be discontinuous, nonconvex, multimodal, noisy-to-evaluate, and so forth ³, especially in control problems ¹⁵. Furthermore, in multiobjective optimization problems, evolutionary algorithms evolve a whole set of solutions that constitutes an estimate of the nondominated set in a single run. In contrast, most of the other algorithms must perform a complete run for finding each point over the set.

In this article, a stochastic model is used to incorporate an uncertainty element on the biological simulation process, thence, making it more realistic in small communities ⁶. The stochastic SIR (SSIR) intends to image this behavior by introducing random “noise” with zero mean into the deterministic system equations ²¹. SSIR also can be used to cope with vaccination, with the intention to the global eradication of a disease, which is, for sure, desirable for the World Health Organization. Vaccination in stochastic models is also considered in several works. For example, in paper ¹³ epidemic models with continuous vaccination strategies are investigated, allowing random fluctuation around the endemic equilibrium, and the transmission rate is analyzed. A mathematical study about optimal vaccination considering a random variable describing the length of the infectious period of a typical infectious is performed in paper ⁴. Finally, epidemic stochastic models including a vaccination variable are considered, for example, in paper ²⁴, which studies the existence, uniqueness, and positivity of the solution and the stability of disease free equilibrium.

The SSIR model is extended here to analyze the effects of impulsive vaccination on the population, which means that in certain time steps, a given percentage of the susceptible people gets vaccinated and becomes part of the recovered population, and the dynamic system keeps autonomous in the time intervals between the consecutive impulse actions of control. It proposes a computer simulation and optimization model: to use the SSIR model with fixed parameters as the dynamic system in a multiobjective impulsive control scheme based on an open-loop continuous-variable dynamic optimization procedure, regarding, simultaneously, with two main targets which conflict with each other: making the probability of eradication as small as possible and making the vaccination campaign as cheap and handy as possible.

The proposed strategy is to be made in two phases. Phase I corresponds to the transient regime control, taking into account the necessity to reduce the number of infected individuals to an acceptable level in a finite time, starting from an equilibrium, for example, considering as decision variables sequences of ratios of susceptible to be vaccinated. And Phase II, corresponding to the permanent regime control, that will act since the result of a solution in Phase I is chosen, to avoid an/another outbreak, considering fixed values of susceptible that will be vaccinated and the time interval between each campaign. The final solution is a set of nondominated points, where one objective can not be reduced without increasing the other one.

The article is organized as follow: Section 2 presents the deterministic SIR model; Section 3 presents the stochastic counterpart model and Section 4 explain formulations on how to adapt SSIR to consider impulsive vaccination via multiobjective optimization. Section 5 presents the numerical simulation results and Section 6 presents final discussions and the conclusions of the work.

2 THE DETERMINISTIC SIR MODEL

The SIR epidemiological model, proposed by Kermack and McKendrick in 1927 ¹⁹, is widely considered to represent infectious diseases ¹⁷. This model examines the spread of a disease in a population during a time-stage, in a average way, classifying the individuals in three states:

• Susceptible(S) - individuals which are not infected, but may become by interaction with infected ones;

• Infectious (I) - individuals which are sick and can transmit the disease to susceptible ones;

• Recovered (R) - individuals which had been sick or received vaccine and are now immune to the contagion.

These states are related by a system of nonlinear differential equations, shown in Equation (2.1). Individuals when born are on the susceptible state, and for every individual that dies, a new susceptible born. Susceptible individuals can be infected by contact transmission, and an infected individual may be cured, becoming recovered and immune.

$\left\{\begin{array}{c}\frac{dS}{dt}=\mu N-\mu S-\frac{\beta SI}{N},S_o\ge 0,\\ \\ \frac{dI}{dt}=\frac{\beta SI}{N}-\gamma I-\mu I,I_o\ge 0,\\ \\ \frac{dR}{dt}=\gamma I-\mu R,R_o\ge 0,\end{array}\right.$ (2.1)

in which N is the number of individuals, $N=S\left(t\right)+I\left(t\right)+R\left(t\right)$, µ is the rate of new susceptible, β is the disease transmission rate from infected to susceptible individuals, and γ is the recovery rate of infected individuals. The set of initial conditions in $t=0$ must be previously known, if the explicit simulation of the system is to be performed. These equations can also be considered in percentages, with $s:S/N,i:I/N,r:R/N$ and the variable r must not be considered.

The basic reproductive number R ₀ is defined in Equation (2.2), which represents the average number of secondaries infections produced by an infected individual.

$R_0=\frac{\beta }{\mu +\gamma }$ (2.2)

The article ¹⁷ shows that the SIR system has an asymptotically stable endemic $(i\ge 0)$ equilibrium iff $R_0>1$, or $\beta \ge \mu +\gamma$, given in Equation (2.3).

$\left(s^{*},i^{*}\right)=\left(\frac{1}{R_0},\frac{\mu }{\beta }(R_0-1)\right)$ (2.3)

It is worth noticing that the percentage equilibrium does not depend on the population size N, and, for extreme parameter values the stable point can describe a situation in which the disease is almost eradicated (if $\beta <<\gamma$ and µ) or in which almost the whole population gets infected (if $\beta >>\gamma$ and µ).

The parameters used in this work are based on the measles epidemic in 1950-68 in the U.K. ², according Table 1. Notice that the stable point is, in percentage, $(s*,i*)\approx (0.0589,0.2403)$, which means that the amount of infected individuals in the equilibrium is very significantly endemic, and it should be necessary to consider policies to control the disease. In fact, the basic reproductive number for this case is $R_0\approx 17$. For this reason, a control action is necessary to reduce the total number of infected individuals.

Table 1:
Parameters of SIR Model.

	µ	β	γ
Values	0.01	0.98	0.04

3 THE STOCHASTIC SIR MODEL

The development of an epidemic in nature underlies many influences which superpose each other and which are therefore practically not predictable. The stochastic SIR (SSIR) model is used intending to better image the behavior of Susceptible, Infected and Recovered by introducing random “noise” into the deterministic system. It is assumed that the epidemic processes, such as infection or death due to infection, are deterministic, but every other process controlling the demography are stochastic. The weapon of choice is usually, of course, the Brownian motion, since it has mean zero, it is symmetric regarding negative and positive deviations. It reflects and describes almost every random situation that can be observed in nature or socio-economics ^{5), (20}.

The stochastic model is, at first, a continuous time Markov process $\left\{\left(S_t,I_t,R_t\right):t\in [0,T_f]\right\}$, with values in ${\mathrm{\mathbb{Z}}}_{+}^3$. The rates from above become the conditional transition rates of the stochastic process (S, I, R), for example: $P(S_{t+\Delta t}=s+1|S_t=s)=\mu N\Delta t+o(\Delta t)$. Each increment of the process is then decomposed into the sum of the expected value of the increment and a sum of independent centered increments among ∆Z _i , relative to each term in the deterministic model:

$\left\{\begin{array}{l}\Delta S=S_{k+1}-S_k=\left(\mu (N-S_k)-\beta \frac{S_k{I}_k}{N}\right)\Delta t+\Delta Z_1-\Delta Z_2,\\ \Delta I=I_{k+1}-I_k=\left(\beta \frac{SI}{N}-(\gamma +\mu )I_k\right)\Delta t+\Delta Z_2-\Delta Z_3.\end{array}\right.$ (3.1)

Each increment ∆Z _i has zero mean (since the expected population size shall be constant) and its variance is expressed in Equation (3.2). The variance ${\sigma }_1^2$ is made by adding up the variances of births and deaths, and so forth.

$\left\{\begin{array}{c}{\sigma }_1^2={\sigma }^2\left(\Delta Z_1\right)=\mu (N+S)\Delta t,\\ \\ {\sigma }_2^2={\sigma }^2\left(\Delta Z_2\right)=\beta SI/N\Delta t,\\ \\ {\sigma }_3^2={\sigma }^2\left(\Delta Z_3\right)=(\gamma +\mu )I\Delta t.\\ \end{array}\right.$ (3.2)

It is a well known fact that for small time steps ∆t and large population size N, a process with centered Poisson increments is well approximated by the process with Gaussian increments having the same conditional variances ²¹. Therefore, if the dimensionless variables are introduced and each Poisson variables ∆Z _i in Equation 3.1 is replaced by increments of Brownian motion or - mathematically correct - Wiener process dW _i , the Equation (3.1) becomes as the model described in Equation (3.3), which is the SSIR system used in this paper.

$\left\{\begin{array}{c}dS=\left(\mu (N-S)-\beta \frac{SI}{N}\right)dt+G_{1}d{W}_1-G_{2}d{W}_2,S_o\ge 0;\\ \\ dI=\left(\beta \frac{SI}{N}-(\gamma +\mu )I_k\right)dt+G_{2}d{W}_2-G_{3}d{W}_3,I_o\ge 0.\\ \end{array}\right.$ (3.3)

The standard deviations G _i are designed by Equation (3.4).

$G_1=\sqrt{\mu (N+S)},G_2=\sqrt{\beta SI/N},G_3=\sqrt{(\gamma +\mu )I}.$ (3.4)

Therefore, it is clear that, generally, for large values of N, the stochastic model will tend to the deterministic model. In general, the smaller the value of N, the bigger will be the influence of the noise. On the other hand, with increasing values of N, the influence of the stochastic disturbances decreases and the stochastic model converges to the deterministic model, which means that the characteristics and solution of the stochastic system will be the same as those of the deterministic system, if N is sufficiently large. This is also shown by the pictures below.

Each realization of this model is numerically obtained by using the Euler-Maruyama method, the stochastic version of the well known Euler method for solving stochastic differential equations ^{18), (20}, with 10, 000 trials. The parameters used in both simulations are shown in Table 1, and $s_o=0.98,i_o=0.02,N=200$.

This model can then be used to estimate probabilities, which distinguish it from its deterministic counterpart. An important feature in stochastic epidemic models is the possibility of computing the probability of a disease be eradicated.

The probability of eradication (POE) can be simply estimated, via Monte Carlo Simulations ²³, given a final time T _f and a number of execution realizations N _R . Numerically, the POE can be calculated by the ratio between the number of runs whose the final infected number is zero or nearby zero and the number of total runs, according to Equation (3.5).

$\text{POE}=\frac{{\displaystyle \#\{I(T_f)\le ϵ\}}}{{\displaystyle N_R}}$ (3.5)

It is noticeable that for small populations, the disease can be eradicated by chance, which means that the POE is not null. In fact, the population size has influence on the behavior of our system, which is not the case for the deterministic model where the size of the population does not influence the solution. In fish farms or cattle herds for example, a population size could be chosen which maximizes the POE given a certain vaccination factor. If needed, the cattle or fish population would be broken down into partitions of smaller sizes given by the information taken from the picture.

Especially, for small N the effect of the noise is great which can lead to some strange behavior of the system. In general, the number of infectious subjects is typically small at the beginning and for that reason the probability to infect a large percentage of the susceptible part of the population is also small. With the influence of chance it may happen that from a certain time no further susceptibles become infected and the disease is eradicated without any external factors at work. Therefore, just for small population communities it may not be necessary to perform control. By the other side, for high population size, control actions may be necessary when the deterministic equilibrium indicates that it is hard to eradicate the disease by chance.

4 OPTIMIZATION MODEL

It is well known the positive impact of vaccination on the public health and its major effect on the mortality reduction and the population growth around the world ²². In fact, since the nineteenth century, it has been saw a sequence of successful campaigns, for example, to reduce and eradicate hepatitis, yellow fever, polio and measles. Although vaccination against a disease is nowadays one of the most usual techniques to control an epidemic in human people, the organization of the campaigns is still a big challenge to governs, specially in non-developed countries, aimed to eradicate the disease, or keeping it in a reasonable level with a minimal cost.

This work is based on control by pulse vaccination, which means that in certain time steps, a given percentage of the susceptible individuals get vaccinated and becomes part of the recovered population, according to references ^{1), (7), (9), (10), (14}. This approach allows different sizes of pulse control actions at arbitrary time instants. It is a multiobjective control since the target is to optimize both vaccination costs and the spending with infected.

Previous sections have shown that, free of vaccination influences, the deterministic and stochastic SIR systems converge to a equilibrium point which depends on the choice of the parameters (Equation 2.3). If an impulsive vaccination is applied, the formerly continuous system becomes non-continuous with the functions S and R becoming step functions according to the vaccination and the function of I reacting according to the system of impulsive differential equations.

4.1 Background in Multiobjective Optimization

Multiobjective problems mean that there are r-scalar objective-functions $F:{\mathrm{ℝ}}^n\to {\mathrm{ℝ}}^r$, with the intention to be simultaneously minimized, in which r is the total number of objectives of the problem. The solutions of a multiobjective problem are inside the Pareto-optimal set or non-dominated set. A solution is said to be non-dominated if another point that is better for one of the functions, is worse for at least one other function, which means that a solution is only dominated by another point if that second one is better in every considered objective-functions.

Mathematically, in a minimization problem, if y ∈ ℝⁿ denotes the vector of decision variables of the problem and $\mathcal{D}\subset {\mathrm{ℝ}}^n$ the set of feasible solutions y, the Pareto-optimal set $\mathcal{P}\subset \mathcal{D}$ is characterized by:

$\mathcal{P}\triangleq \left\{\overline{y}\in \mathcal{D}|\nexists y\in \mathcal{D}:y\prec \overline{y}\right\}.$ (4.1)

This relation of domination, denoted by $y\prec \overline{y}$, means that given two different decision variable vectors, y and $\overline{y}$ (which lead respectively to the cost vectors F(y) and $F\left(\overline{y}\right)$), then it has $F_i\left(y\right)\le F_i\left(\overline{y}\right),\forall i=1,2,\dots ,r$, and $F_i\left(y\right)<F_i\left(\overline{y}\right)$, for some i.

Evolutionary multiobjective optimization (EMO) algorithms have become an important computational tool for solving problems with multiple conflicting objective functions. This has occurred, in part, because while the evolutionary algorithms evolve a whole set of solutions that constitutes an estimate of the problem Pareto-set in a single run, most of the other algorithms which could be used for performing the same task must perform a complete run for finding each point over the Pareto-set. EMO techniques can deal with objective functions which can be discontinuous, nonconvex, multimodal or noisy-to-evaluate and so forth. These algorithms are heuristics; in other words, it is designed to find feasible solutions to the problem but there is no guarantee for how close those solutions are to the true optimal solution. More extensive information can be found, for example, in paper ⁸.

Genetic algorithms intend to copy the behavior of nature who optimizes using the principles of selection, mutation and recombination, in which a set of solutions is considered in each iteration. The set of these solutions undergoes a selection process, when bad solutions “die” and are not considered anymore, better solutions “survive”, making up a set of parents “with good genes” for the next rounds. From this set of considerably good solutions a new population is created by employing the principles of mutation and recombination. This new set then undergoes the selection process again. Solutions are always compared regarding their state of being dominated by others or being dominant for other solutions. At the end of the process the best solutions, those who are not dominated by any other solution, form the Pareto-front ¹⁶. This paper deals with the question of finding an estimate of the Pareto-optimal set of the problems using a parallel version of the standard NSGA-II (The Non-dominated Sorting Genetic Algorithm) ¹², which is a fast and elitist multiobjective optimization engine, with selection considering both the non-dominance ranking and the crowding distance.

4.2 Applying pulse vaccination

Depending on the set of given parameters (µ, β, γ) on a SIR system (deterministic or stochastic), impulsive vaccinations can be applied if the control of the disease is wished, which means that the campaigns, for example, should keep the percentage of infected people under a given level. Other than this, it is often desired to reduce the significant peak of I, the maximum of the number of infected people within the population.

Applying vaccination campaigns mean that in certain time steps a given percentage of the susceptible population gets vaccinated against the disease and becomes part of the recovered population. Then, the formerly continuous systems become non-continuous with the functions S and R becoming step functions according to the vaccination and the function of I and reacting according to the system of differential equations.

To be precise, according to ⁷, it means that the time in the closed interval $[0,T_f]$, being T _f the final time, has to be partitioned in a set of instants previously fixed $\Gamma =\{{\tau }_0,...,{\tau }_M\}$, such that: ${\tau }_0=0,{\tau }_M=T_f$ and ${\tau }_{k+1}-{\tau }_k=\Delta T$. The time between each τ _k and τ _k+1 is considered as a stage of the problem. Consider the percentage of the susceptible individuals that will be vaccinated during the campaigns: $P=\left\{p\left[0\right],...,p[M-1]\right\}$, such that $p\left[k\right]=p\left({\tau }_k\right)$, for each τ _k in Γ.

Since the vaccinations act impulsively only on the ratio of the susceptible people, and considering the time ${\tau }_k^{+}$ as the time instant “just before” the impulsive action in τ _k , thus:

$\left\{\begin{array}{l}s\left({\tau }_k^{+}\right)=s\left({\tau }_k\right)(1-p[k\left]\right),\\ i\left({\tau }_k^{+}\right)=i\left({\tau }_k\right),\\ k=0,1,...,M-1.\end{array}\right.$ (4.2)

Each Initial Value Problem (ref. Equations (2.1) or (3.3)) is valid during the time in which there is no control action (no vaccination), the system presents its autonomous dynamics within such time intervals, and a new initial condition is established in each optimization stage according to Equation (4.2), linking each stage to the next one. The set of initial conditions in $t=0$ must be previously known, if the explicit simulation of the system is to be performed. If convenient, the percentage of the susceptible that will be vaccinated during the campaigns may be fixed: $p\left[k\right]=p$, for each τ _k in $k=0,1,...,M-1$.

4.3 Multiobjective formulation

It is clear that higher vaccination ratios will in general lead to a smaller amount of infectious persons within the population and may even lead to extinction of all infected individuals. Meanwhile, different vaccination campaigns may lead to the same probability of eradication. If the intention is keeping the percentage of infected people under a given level, or to reduce the significant peak of the number of infected people within the population, it is desirable to be made trying to spend as lowest as possible in vaccination campaigns. This can be achieved, for example, by vaccinating much more people at the beginning of the campaign while later it is only needed to vaccinate few people to further control the disease and keep it under a certain level.

Multiobjective optimization offers a whole set of optimal solutions regarding the total number of infected people as well as the costs of the vaccination campaign. Responsible managers can then decide for one of these optimal campaigns, according to other variables like budget, media influence, society view of the diseases, guidelines by law or politics and so on. The diversity of the Pareto-optimal solutions will permit to the public manager choose the policy that can better be applied to each present situation, with or without a greater prevention cost, but always spending as less as possible in each case.

Here, similar to paper ¹⁰, the optimization is made in two phases:

• The first phase corresponds to the transient regime control, taking into account the necessity to reduce the number of infected individuals to an acceptable level in a finite time, starting from an equilibrium, for example, considering as decision variables sequences of ratios of susceptible to be vaccinated.

• An unique solution of Phase I will be selected and implemented, and Phase II, corresponding to the permanent regime control, will act since the result of this solution, considering as decision variables fixed value of susceptible that will be vaccinated and the time interval between each campaign, to avoid an/another outbreak.

The intention is to find a set of equally-optimal solutions, called the Pareto-optimal set, or the set of non-dominated points.

Since a stochastic dynamic model is considered here, the application of a same vaccination campaign may have different results. For example, the same vaccination campaign may lead to extinction of the disease during the considered time horizon while in other runs, the disease may not be eradicated and survive. If a large number of simulations of the same kind is run and the ratio of number of outcomes with eradication contra the number of outcomes without eradication considered, the probability of eradication (POE) for a certain vaccination campaign can be approximated, which will be an objective func∀tion. The other objective is to minimize the cost of the control policy.

4.3.1 Formulation of the permanent regime control

The optimization variables of the problem are:

• the time interval between the campaigns: ∆T,

• the fixed percentage of the susceptible that will be vaccinated during the campaigns: $p=p\left[k\right]$ , for $k=0,1,...,M-1$ .

The optimal vaccination campaign planning is formulated here as a multiobjective optimization problem with two objective functions that should be minimized:

• F₁ - the probability of eradication during the optimization horizon,

• F₂ - the total mean cost with the campaign (sum of vaccinated percentages).

The model is subject to the follow constraints:

• $0<\Delta T<T_f$ ,

• $p_{min}<p<p_{max}$ ,

for given percentage ratios p _min and p _max .

The optimization problem is then showed in Equations (4.3) and (4.4):

$\underset{p,\Delta T}{\min }\left\{\begin{array}{l}{F}_1=\frac{{\displaystyle \#\{I(T_f)\le ϵ\}}}{{\displaystyle N_R}}\\ F_2=\sum _{k=0}^{M-1}p\left[k\right]\end{array}\right.$ (4.3)

subject to:

$\left\{\begin{array}{l}\frac{{\displaystyle dS}}{{\displaystyle dt}}=\mu N-\mu S-\frac{{\displaystyle \beta IS}}{{\displaystyle N}},S_o\ge 0;\\ \frac{{\displaystyle dI}}{{\displaystyle dt}}=\frac{{\displaystyle \beta IS}}{{\displaystyle N}}-\gamma I-\mu I,I_o\ge 0;\\ s\left({\tau }_k^{+}\right)=s\left({\tau }_k\right)(1-p);\\ i\left({\tau }_k^{+}\right)=i\left({\tau }_k\right);\\ k=0,1,...,M-1;\\ \begin{array}{l}0<\Delta T<T_f\\ p_{min}<p<p_{max}\end{array}\end{array}\right.$ (4.4)

4.3.2 Formulation of the transient regime control

The parameters of the optimization model are:

• the number of vaccination campaigns (stages): M,

• the fixed vector of time instants between the campaigns: $\Gamma =\{{\tau }_0,\dots ,{\tau }_M\}$ .

The optimization variables of the problem are:

• the percentages of the susceptible individuals that will be vaccinated during the campaigns: $P=\left\{p\right[1],\dots ,p[M\left]\right\}$ , such that $p\left[k\right]=p\left({\tau }_k\right)$ , for each τ _k in Γ.

The constraints are modeled as:

• Each vaccination ratio p[k] has to follow the rule: $0<p_{min}\le p\left[k\right]\le p_{max}\le 1.0$ ,

for given percentage ratios p _min and p _max .

The vaccination campaign planning is formulated as a multiobjective optimization problem with two objective functions that should be minimized:

• F₁ - the probability of eradication during the optimization horizon,

• F₂ - the total mean cost with the campaign (sum of vaccinated percentages).

The optimization problem is then showed in Equations (4.5) and (4.6):

$\underset{P}{\min }\left\{\begin{array}{l}{F}_1=\frac{{\displaystyle \#\{I(T_f)\le ϵ\}}}{{\displaystyle N_R}}\\ F_2=\sum _{k=0}^{M-1}p\left[k\right]\end{array}\right.$ (4.5)

subject to:

$$\begin{gathered} \left\{\begin{array}{l}\frac{{\displaystyle dS}}{{\displaystyle dt}}=\mu N-\mu S-\frac{{\displaystyle \beta IS}}{{\displaystyle N}},S\left(0\right)=S_o\ge 0;\\ \frac{{\displaystyle dI}}{{\displaystyle dt}}=\frac{{\displaystyle \beta IS}}{{\displaystyle N}}-\gamma I-\mu I,I\left(0\right)=I_o\ge 0;\\ t\in ({\tau }_k^{+},{\tau }_{k+1}];\\ s\left({\tau }_k^{+}\right)=s\left[k^{+}\right]=s\left[k\right](1-p[k\left]\right); \\ i\left({\tau }_k^{+}\right)=i\left[k^{+}\right]=i\left[k\right]; \\ k=0,1,\dots ,M-1;\\ 0<p_{min}\le p\left[k\right]\le p_{max}\le 1.0;\end{array}\right. \end{gathered}$$ (4.6)

5 NUMERICAL RESULTS

This paper deals with the statement of finding an estimate of the Pareto-optimal set of both permanent and transient regime control problem. The parameters considered here are: $T_f=360u.t.,p_{min}=0.4,p_{max}=0.8$. The SIR model parameters values are in Table 1.

First of all, for better comparison, Figures (1a) and (1b) show the evolution of the population of susceptible S, infectious I and recovered R considering a few simulation runs of the stochastic model and the evolution came from the deterministic SIR model, respectively, without control. On the one hand, notice that even though the behavior of the epidemic is never the same, the stochastic SIR seems to converge, in an average way, to the same stable point of its deterministic counterpart. Furthermore, we can see that, according to Section 2, the stable point without vaccination is, in percentage, $(s^{*},i^{*})\approx (0.0589,0.2403)$. It means that the amount of infected individuals in the equilibrium is very significantly endemic (around 24% of the population). In addition, there is a disease peak of approximately 80% of the population infected in the very beginning of the evolution time. For these reasons, a control action is necessary to reduce the total number of infected individuals.

Figure 1:
Evolution using the stochastic and deterministic SIR without vaccination.

Therefore, a parallel version of the standard NSGA-II is considered, distributing the evaluation of the objective function. Each algorithm is tested in 30 independent executions, and each SSIR is simulated in $N_R=150$ trial runs in 8 cores of a standard workspace. The population size is $\mathit{\text{popsize}}=40$, the number of generations is $nger=100$, the probability of crossover is $p_c=0.8$, the probability of mutation is $p_m=0.05$.

5.1 Numerical results for the permanent regime control

The result come from the multiobjective optimization procedure is presented in Figures 2 and 3. They show the nondominated set in the objective space and in the parameter space as well. Doing the analysis of this figure, the responsible will be able to perform a more confidant decision over the Pareto-optimal set, by using the relation between the POE and the cost of the vaccination campaigns.

Figure 2:
Nondominated fronts (costs with vaccination as a function of the probability of eradication) for the permanent regime control.

Figure 3:
Nondominated variables (fixed percentage of vaccinated people and time interval between the campaigns) as a function of their probability of eradication came from the multiobjective optimization for the permanent regime control.

Figure 2a shows the combined Pareto front considering the probability of eradication versus the cost with vaccination came from the 30 executions of the NSGA-II and Figure 2b shows the final nondominated set obtained after applying a nondominance procedure. Note that the leftmost campaigns within the Pareto front have POE equal to zero, according to the fact that their integral of infected persons is highest. The same holds for the rightmost campaigns from the Pareto front: since all of them share a POE equal to one. In fact, these two graphics are similar, which means that the optimization procedure is robust, once all the executions converge approximately to the same final result. Also note that the Pareto set is composed by two regions: the first one, related to the probabilities of eradication between 0 and 0.4, has greater slope than the other one, related to the probabilities of eradication between 0.4 and 1. This can be seem as a measure of the big effort necessary to start the process of controlling and eradicating the disease.

The parameter space is shown in Figure 3a, which presents a 3-D graphic of the nondominated vaccination campaigns (the fixed percentage of vaccinated susceptible people and the time interval between the campaigns) as a function of their probability of eradication. Figure 3b shows the 2-D counterpart of the above considering the fixed percentage of vaccinated people and its probability of eradication. It can be noted the compromise relation between the number of campaigns and the percentage of vaccines that will be applied, when small percentages is related to small probabilities of eradication end large percentages related to large probabilities.

It will now have a closer look at how different strategies influence the population sizes of susceptibles and infectious. Figures 4 and 5 show the result of the application of two distinct control policies referent to different probabilities of eradication. Figures 4a and 4b show the result referent to the probability of eradication equal to 0.3 and Figures 5a and 5b show the result referent to the probability of eradication equal to 0.8. It shows the optimal campaign (the fixed percentage of vaccinated people and time interval between the campaigns) and the system evolution (S(t), I(t) and R(t)) for these two distinct strategies. Note that according to the percentage of vaccinated persons the number of susceptibles changes, giving a non-continuity in the functions describing the susceptibles and recovered. Other than that, the infected part of the population reacts according to the laws of the model.

Figure 4:
Simulation considering two nondominated solutions related to probability of eradication 0.3 in multiobjective optimization for the permanent regime control.

Figure 5:
Simulation considering two nondominated solutions related to probability of eradication 0.8 in multiobjective optimization for the permanent regime control.

Finally, it can be noticeable that the optimization is influencing the cost minimization. In fact, the graphics show that higher vaccination ratios will, in general (the simulation depends on randomness) lead to a smaller amount of infectious people within the population and may finally even lead to extinction of all infected persons. Comparing with the autonomous case, with no vaccination, when around 24% of the population is infected, the graphs show that, in the case of less control, with the vaccination rate of 30%, the percentage of infected is below 10% of the population. On the other hand, in the case of greater control, with a vaccination rate of around 80%, the percentage of vaccinated persons is very close to zero, which corresponds to silencing the disease. Additionally, notice that the peak of the disease has practically disappeared with the proposed control actions.

5.2 Numerical results for the transient regime control

The result come from the multiobjective optimization procedure is presented in Figure 6 and 7. They show the nondominated set in the objective space and in the parameter space as well. Figure 6a shows the combined Pareto front considering the probability of eradication versus the cost with vaccination came from the 30 executions of the NSGA-II and Figure 6b shows the final nondominated set obtained after applying a nondominance procedure. Note that these two graphics are similar, but the convergence is slow, once the number of optimization variables are an order of magnitude higher. Note also that they are similar to the nondominated set came from the optimal permanent regime control.

Figure 6:
Nondominated fronts (costs with vaccination as a function of the probability of eradication) for the transient regime control.

Figure 7:
Nondominated variables (campaigns consisting in a percentage of vaccinated people for each time interval) as a function of their probability of eradication came from the multiobjective optimization for the transient regime control.

The parameter space is shown in Figure 7a, which presents a 3-D graphic of the nondominated vaccination policies, consisting in, for each nondominated individual (shown as 0 to 100), the set of optimal percentage of vaccinated people sequence for each time interval. Figure 7b shows a graphic with the nondominated vaccination policies (shown as 0 to 100) as a function of their probability of eradication.

Figures 8 and 9 show the result of the application of two distinct control policies referent to different probabilities of eradication. Figures 8a and 8b show the result referent to the probability of eradication equal to 0.3 and Figures 9a and 9b show the result referent to the probability of eradication equal to 0.8. It shows the optimal campaign (the sequence of the percentage of vaccinated people and time interval between the campaigns) and the system evolution (S(t), I(t) and R(t)) for these two distinct policies.

Figure 8:
Simulation considering two nondominated solutions related to probability of eradication 0.3 in multiobjective optimization for the transient regime control.

Figure 9:
Simulation considering two nondominated solutions related to probability of eradication 0.8 in multiobjective optimization for the transient regime control.

Also here, it can be noticeable that the optimization is influencing the cost minimization. The difference is that, as in many practical cases, the rate of vaccination here decreases over the time specified for control. This implies a time-evolution of the system in which two distinct behaviors are observed: a drastic reduction of the infected people in the phase in which the vaccination is applied and a gradual increase in the number of infected people in the phase in which the vaccination rate reduces or there is no more control. In fact, the number of infected people returns to the equilibrium in both cases, with around 24% of the population is infected applying vaccination policies with a probability of eradication 80% and 30%. In any case, also notice that the peak of the disease has practically disappeared with the proposed control actions.

Finally, it should be pointed that depending on the parameters chosen for the model (like population size and time horizon) for NSGA-II (like population size in each generation, number of generations, and number of decision variables), the time for optimization is huge, even implementing parallelization. Just one function evaluation takes several minutes because for calculating the POE, many realizations of the stochastic model have to be run for computing the ratio of runs with eradication compared to the total number of runs. Because every individual in the NSGA-II has to undergo this function evaluation in each generation, that time has to be multiplied by the number of individuals and the number of generations that the NSGA-II uses. Sorting of the individuals and computing of the crowding distance which is needed for that also takes some time resulting in several days which are required for one entire optimization run.

6 CONCLUSIONS

This work considered a quantitative study of a general epidemic control via vaccination. It proposed a framework to analyze the effects of impulsive vaccinations to epidemics control in two phases: the transient regime, reducing the number of infected individuals to an acceptable level in a finite time, and the permanent regime, acting with fixed vaccinations to avoid another outbreak. The stochastic and deterministic SIR model is considered with parameters chosen to obtain a very significantly endemic infected individuals in the equilibrium, based on the measles epidemic in 1950-68 in the U.K. data. A multiobjective optimization problem was proposed, regarding, simultaneously, with two main targets which conflict with each other: making the probability of eradication as small as possible and making the vaccination campaign as cheap and handy as possible. An estimate of the Pareto-optimal set of both problems was found using a parallel version of the NSGA-II.

The proposed scheme gives back the sequence of people that should be vaccinated The final result shows nondominated policies that can be a useful guidance for public managers to decide which is the best procedure to be taken depending on the present situation.

Numerical results showed that the control action had been needed to reduce the total number of infected individuals, and the proposed model and methodology do contribute to this reduction. This can be seen as a measure of the significant effort necessary to start controlling and eradicating the disease. n the permanent regime control, the graphs showed that, in the case of less vaccination control, the percentage of infected was below 10% of the population and in the case of greater vaccination control, the percentage of vaccinated persons was very close to zero, which corresponds to silencing the disease. In the transient regime control, the vaccination rate decreased over time, leading to a time-evolution with two distinct behavior: a drastic reduction of the infected people when the vaccination is being applied and a gradual increase in the number of infected people when the vaccination rate reduces. Additionally, in both cases, the peak of the disease has practically disappeared with the proposed control actions.

Future works may consider real data case studies and other features can be considered, such as isolation or social distancing and the efficiency of vaccines.

Another segment would be to investigate the optimization machinery to reduce the computational cost of the execution, not lose the diversity of the population and the quality of the set of non-dominated solutions.

This study can also suggest owners of large cattle herds or fish farms as to which population size to choose in case of the outbreak of an epidemic.

Acknowledgements

The authors would thank to the Brazilian Agencies CNPq grant 423864/2016-5, FAPEMIG grant APQ-02778-17, and to CEFET-MG and IAESTE for financial support.

REFERENCES

• ¹
  K. Adam & R.T.N. Cardoso. Optimal Multiobjective Pulse Vaccination Campaigns in Stochastic SIR Model. In “Congresso Nacional de Matemática Aplicada e Computacional”, volume 1. SBMAC (2019).
• ²
  R.M. Anderson & R.M. May. “Infectious diseases of humans: dynamics and control”. Oxford University Press (1992).
• ³
  T. Bäck. “Evolutionary algorithms in theory and practice: evolution strategies, evolutionary programming, genetic algorithms”. Oxford University Press, USA (1996).
• ⁴
  F. Ball & O. Lyne. Optimal vaccination policies for stochastic epidemics among a population of households. Mathematical Biosciences, 177 (2002), 333-354.
• ⁵
  F. Brauer P., Driessche & J. Wu. “Mathematical Epidemiology (Lecture Notes in Mathematics)”. Springer (2008).
• ⁶
  T. Britton. Stochastic epidemic models: a survey. Mathematical biosciences, 225(1) (2010), 24-35.
• ⁷
  R. Cardoso & R. Takahashi. Solving impulsive control problems by discrete-time dynamic optimization methods. Tendências em Matemática Aplicada e Computacional, 9(1) (2008), 21-30.
• ⁸
  C.A.C. Coello G.B. Lamont & D.A.V. Veldhuizen. “Evolutionary Algorithms for Solving Multi-Objective Problems”. Springer (2007).
• ⁹
  A. da Cruz R. Cardoso & R. Takahashi. Multiobjective dynamic optimization of vaccination campaigns using convex quadratic approximation local search. In “Evolutionary Multi-Criterion Optimization”. Springer (2011), pp. 404-417.
• ¹⁰
  A. da Cruz R. Cardoso & R. Takahashi. Multiobjective synthesis of robust vaccination policies. Applied Soft Computing, 50 (2017), 34-47.
• ¹¹
  A.R. da Cruz R.T. Cardoso & R.H. Takahashi. Multi-objective Design with a Stochastic Validation of Vaccination Campaigns. IFAC Proceedings Volumes, 42(2) (2009), 289 - 294. doi:https://doi.org/10. 3182/20090506-3-SF-4003.00053. URL http://www.sciencedirect.com/science/article/ pii/S1474667015367161 14th IFAC Workshop on Control Applications of Optimization.
  » https://doi.org/https://doi.org/10. 3182/20090506-3-SF-4003.00053» http://www.sciencedirect.com/science/article/ pii/S1474667015367161
• ¹²
  K. Deb A. Pratap, S. Agarwal & T. Meyarivan. A fast and elitist multiobjective genetic algorithm: NSGA-II. Evolutionary Computation, IEEE Transactions on, 6(2) (2002), 182-197.
• ¹³
  C. Deng & H. Gao. Stability of SVIR system with random perturbations. International Journal of Biomathematics, 5(04) (2012).
• ¹⁴
  A.C.S. Dusse & R.T.N. Cardoso. Using a Stochastic SIR Model to Design Optimal Vaccination Campaigns via Multiobjective Optimization. In “International Symposium on Mathematical and Computational Biology”, volume 1. BIOMAT (2019).
• ¹⁵
  W.R. Esposito & C.A. Floudas. Deterministic global optimization in nonlinear optimal control problems. Journal of Global Optimization, 17(1-4) (2000), 97-126.
• ¹⁶
  D.E. Goldberg. “Genetic Algorithms in Search, Optimization and Machine Learning”. Addison Wesley Longman, Inc. (1989).
• ¹⁷
  H. Hethcote. The mathematics of infectious diseases. Preprints of the IFAC Workshop on Control Applications of Optimisation, 42(4) (2000), 599-653.
• ¹⁸
  D. Higham. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM review, (2001), 525-546.
• ¹⁹
  W. Kermack & A. McKendrick. A contribution to the mathematical theory of epidemics. Proceedings of the Royal Society of London Series A Mathematical and Physical Sciences, 115 (1927), 700-721.
• ²⁰
  P.E. Kloeden & E. Platen. “Numerical Solution of Stochastic Differential Equations (Applications is Mathematics)”. Springer (1992).
• ²¹
  R. Kuske L. Gordillo & P. Greenwood. Sustained oscillations via coherence resonance in SIR. Journal of theoretical biology, 245(3) (2007), 459-469.
• ²²
  S. Plotkin W. Orenstein & P. Offit. “Vaccines”, volume 707. Elsevier Health Sciences, 5th ed. (2008).
• ²³
  R.Y. Rubinstein & D.P. Kroese. “Simulation and the Monte Carlo method”, volume 707. Wiley-interscience (2008).
• ²⁴
  E. Tornatore S. Buccellato & P. Vetro. On a stochastic disease model with vaccination. Rendiconti del Circolo Matematico di Palermo, 55(2) (2006), 223-240.

Publication Dates

• Publication in this collection
  02 July 2021
• Date of issue
  May-Aug 2021

History

• Received
  28 Dec 2019
• Accepted
  05 Jan 2021