key: cord-241596-vh90s8vi
authors: Libotte, Gustavo Barbosa; Lobato, Fran S'ergio; Platt, Gustavo Mendes; Neto, Antonio Jos'e da Silva
title: Determination of an Optimal Control Strategy for Vaccine Administration in COVID-19 Pandemic Treatment
date: 2020-04-15
journal: nan
DOI: nan
sha: 
doc_id: 241596
cord_uid: vh90s8vi

During decades, mathematical models have been used to predict the behavior of physical and biologic systems, and to define strategies aiming the minimization of the effects regarding different types of diseases. In the present days, the development of mathematical models to simulate the dynamic behavior of novel coronavirus disease (COVID-19) is considered an important theme due to the quantity of infected people worldwide. In this work, the aim is to determine an optimal control strategy for vaccine administration in COVID-19 pandemic treatment considering real data from China. For this purpose, an inverse problem is formulated and solved in order to determine the parameters of the compartmental SIR (Susceptible-Infectious-Recovered) model. To solve such inverse problem, the Differential Evolution (DE) algorithm is employed. After this step, two optimal control problems (mono- and multi-objective) to determine the optimal strategy for vaccine administration in COVID-19 pandemic treatment are proposed. The first consists of minimizing the quantity of infected individuals during the treatment. The second considers minimizing together the quantity of infected individuals and the prescribed vaccine concentration during the treatment, i.e., a multi-objective optimal control problem. The solution of each optimal control problems is obtained using DE and Multi-Objective Differential Evolution (MODE) algorithms, respectively. The results regarding the proposed multi-objective optimal control problem provides a set of evidences from which an optimal strategy for vaccine administration can be chosen, according to a given criterion.

In the last decades, countless mathematical models used to evaluate the spread and control of infectious diseases have been proposed. These models are very important in different fields, such as policy making, emergency planning and risk assessment, definition of control-programs, promoting the improvement of various health-economic aspects (Al-Sheikh, 2013) . In general, such models aim to describe a state of infection (susceptible and infected) and a process of infection (the transition between these states) by using compartmental relations, i.e., the population is divided into compartments by taking assumptions about the nature and time rate of transfer from one compartment to another (Trawicki, 2017; Blackwood and Childs, 2018) . One can cite several studies using models for measles vaccination (Bauch et al., 2009; Widyaningsih et al., 2018) , HIV/AIDS (Mukandavire et al., 2009) , tuberculosis (Bowong and Kurths, 2010) , dengue (Weiss, 2013) , pertussis epidemiology (Pesco et al., 2014) , Alzheimer (Ebrahimighahnavieh et al., 2020) , among others.

Recently, the world has been experiencing the dissemination of a new virus, referred to as COVID-19 (Coronavirus disease 2019). COVID-19 is an infectious disease emerged from China in December 2019, that has rapidly spread around in many other countries worldwide (Gorbalenya et al., 2020; World Health Organization, 2020 (accessed April 8, 2020 . The common symptoms are severe respiratory illness, fever, cough, and myalgia or fatigue, especially at the onset of illness (Huang et al., 2020) . The transmission may happen person-to-person, through direct contact or droplets (Chan et al., 2020; Li et al., 2020; Riou and Althaus, 2020) .

Since the COVID-19 outbreak in Wuhan City in December of 2019, various computational model-based predictions have been proposed and studied. Lin et al. (2020) proposed a Susceptible-Exposed-Infectious-Removed (SEIR) model for the COVID-19 outbreak in Wuhan. These authors considered some essential elements including individual behavioral response, governmental actions, zoonotic transmission and emigration of a large proportion of the population in a short time period. Benvenuto et al. (2020) proposed the Auto Regressive Integrated Moving Average (ARIMA) model to predict the spread, prevalence and incidence of COVID-2019. Roda et al. (2020) used a Susceptible-Infectious-Removed (SIR) model to predict the COVID-19 epidemic in Wuhan after the lockdown and quarantine. In this study, these authors demonstrate that non-identifiability in model calibrations using the confirmedcase data is the main reason for such wide variations. Prem et al. (2020) proposed a SEIR model to simulate the spread of COVID-19 in Wuhan city. In this model, all demographic changes in the population (births, deaths and ageing) were ignored. The simulations showed that control measures aimed at reducing social mixing in the population can be effective in reducing the magnitude and delaying the peak of the COVID-19 outbreak.

In order to evaluate the global stability and equilibrium point of these models, Li and Muldowney (1995) studied a SEIR model with nonlinear incidence rates in epidemiology, in terms of global stability of endemic equilibrium. Al-Sheikh (2013) evaluated a SEIR epidemic model with a limited resource for treating infected people. For this purpose, the existence and stability of disease-free and endemic equilibrium were investigated. Li and Cui (2013) studied a SEIR model with vaccination strategy that incorporates distinct incidence rates for exposed and infected populations. These authors proved the global asymptotical stable results of the disease-free equilibrium. Singh et al. (2017) developed a simple and effective mathematical model for transmission of infectious diseases by taking into consideration the human immunity. This model was evaluated in terms of local stability of both disease free equilibrium and disease endemic equilibrium. Widyaningsih et al. (2018) proposed a SEIR model with immigration and determined the system equilibrium conditions. Kim et al. (2019) developed a Coxian-distributed SEIR model considering an empirical incubation period, and an stability analysis was also performed.

In order to reduce the spread of dissemination of COVID-19 worldwide, various procedures have been adopted. As mentioned by Zhai et al. (2020) and Wei et al. (2020) , quarantine and isolation (social-distancing) can effectively reduce the spread of COVID-19. In addition, wearing masks, washing hands and disinfecting surfaces contribute to reducing the risk of infection. According to the U.S. Food and Drug Administration, there are no specific therapies to COVID-19 treatment. However, treatments including antiviral agents, chloroquine and hydroxychloroquine, corticosteroids, antibodies, convalescent plasma transfusion and radiotherapy are being studied .

As alternative to these treatments, the use of drug administration (vaccine) arises as an interesting alternative to face this pandemic. It must be emphasized that there is currently no vaccine for COVID-19, but there is a huge effort to develop a vaccine in a record time, which justifies the present study (Lurie et al., 2020) . Mathematically, the determination of optimal protocol for vaccine administration characterizes an Optimal Control Problem (OCP). This particular optimization problem consists in the determination of control variable profiles that minimize (or maximize) a given performance index (Bryson and Ho, 1975; Biegler et al., 2002) . In order to solve this problem, several numerical methods have been proposed (Bryson and Ho, 1975; Feehery and Barton, 1996; Lobato, 2004; Lobato et al., 2016) . These methods are classified according to three broad categories: direct optimization methods, Pontryagin's Maximum Principle (PMP) based methods and HJB-based (Hamilton-Jacob-Bellman) methods. The direct approach is the most traditional strategy considered to solve an OCP, due to its simplicity. In this approach, the original problem is transformed into a finite dimensional optimization problem through the parametrization of control or parametrization of control and state variables (Feehery and Barton, 1996) .

From an epidemiological point of view, Neilan and Lenhart (2010) proposed an optimal control problem to determine a vaccination strategy over a specific period of time so as to minimize a cost function. In this work, the propagation of a disease is controlled by a limited number of vaccines, while minimizing a percentage of the overall number of dead people by infection, and a cost associated with vaccination. Biswas et al. (2014) studied different mathematical formulations for an optimal control problem considering a Susceptible-Exposed-Infectious-Recovered model. For this purpose, these authors evaluated the solution of such problems when mixed state control constraints are used to impose upper bounds on the available vaccines at each instant of time. In addition, the possibility of im-posing upper bounds on the number of susceptible individuals with and without limitations on the number of vaccines available were analyzed. The optimal control theory was applied to obtain optimal vaccination schedules and control strategies for the epidemic model of human infectious diseases.

In this word, the objective is to determine an optimal control strategy for vaccine administration in COVID-19 pandemic treatment considering real data from China. In order to determine the parameters that characterizes the proposed mathematical model (based on the compartmental SIR model), an inverse problem is formulated and solved considering the Differential Evolution (DE) algorithm (Storn and Price, 1997; Price et al., 2005) . After this step, two optimal control problems (mono-and multi-objective) used to determine the optimal strategy for vaccine administration in COVID-19 pandemic treatment are proposed. The mono-objective optimal control problem considers minimizing the quantity of infected individuals during the treatment. On the other hand, the multi-objective optimal control problem considers minimizing together the quantity of infected individuals and the prescribed vaccine concentration during the treatment. To solve each problem, DE and Multi-Objective Differential Evolution (MODE) algorithms (Lobato and Steffen, 2011) are employed, respectively.

This work is organized as follows. Section 2 presents the description of the mathematical model considered to represent the evolution of COVID-19 pandemic. In Section 3, the general aspects regarding the formulation and solution of an OCP is presented. A brief review on DE and its extension to deal with multi-criteria optimization is presented in Section 4. In Section 5, the proposed methodology is presented and discussed. The results obtained using such methodology are presented in Section 6. Finally, the conclusions are outlined in Section 7.

In the specialized literature, various compartmental models used to represent the evolution of an epidemic can be found (Forgoston and Schwartz, 2013; Pesco et al., 2014; Shaman et al., 2014; Cooper et al., 2016; Azam et al., 2020) . The study of these models is very important to understand an epidemic spreading mechanism and, consequently, to investigate the transmission dynamics in population (Forgoston and Schwartz, 2013) . As mentioned by Keeling and Rohani (2007) , these compartmental models can be divided into two groups: i) population-based models and ii) agentbased or individual-based models. In turn, the first one can be subdivided into deterministic or stochastic (considering continuous time, ordinary differential equations, partial differential equations, delay differential equations or integrodifferential equations) or discrete time (represented by difference equations). The second class can be subdivided into usually stochastic and usually discrete time.

In the context of population-based models, the deterministic modeling can be represented, in general, by the interaction among susceptible (denoted by S -an individual which is not infected by the disease pathogen), exposed (denoted by E -an individual in the incubation period after being infected by the disease pathogen, and with no visible clinical signs), infected/infectious (denoted by E -an individual that can infect others) and recovered individuals (denoted by R -an individual who survived after being infected but is no longer infectious and has developed a natural immunity to the disease pathogen). Considering a population of size N, and based on the disease nature and on the spreading pattern, the compartmental models can be represented as (Keeling and Rohani, 2007; Hethcote, 2000) :

• Susceptible-Infected (SI): population described by groups of Susceptible and Infected;

• Susceptible-Infected-Recovered (SIR): population described by groups of Susceptible, Infected and Recovered;

• Susceptible-Infectious-Susceptible (SIS): population also described by groups of Susceptible and Infected. In this particular case, recovering from some pathologies do not guarantee lasting immunity. Thus, individuals may become susceptible again;

• Susceptible-Exposed-Infectious-Recovered (SEIR): population described by groups of Susceptible Exposed, Infected and Recovered.

It is important to mention that in all these models, terms associated with birth, mortality and vaccination rate can be added. In addition, according to Keeling and Rohani (2007) and Hethcote (2000) , these models can include: i) time-dependent parameters to represent the effects of seasonality; ii) additional compartments to model vaccinated and asymptomatic individuals, and different stages of disease progression; iii) multiple groups to model heterogeneity, age, spatial structure or host species; iv) human demographics parameters, for diseases where the time frame of the disease dynamics is comparable to that of human demographics. Human demographics can be modeled by adopting constant immigration rate, constant per capita birth and death rates, density-dependent death rate or disease-induced death rate. Thus, the final model is dependent on assumptions taken during the formulation of the problem.

In this work, the SIR model is adopted, in order to describe the dynamic behavior of during the COVID-19 epidemic in China. The choice of this model is due to the study conducted by Roda et al. (2020) . These authors demonstrated that the SIR model performs more adequately than the SEIR model in representing the information related to confirmed case data. For this reason, the SIR model will be adopted here. The schematic representation of this model is presented in Fig. 1 . Mathematically, this model has the following characteristics:

• An individual is susceptible to an infection and the disease can be transmitted from any infected individual to any susceptible individual. Each susceptible individual is given by the following relation:

where t is the time, β and µ represents the probability of transmission by contact and per capita removal rate, respectively. In turn, S 0 is the initial condition for the susceptible population.

• Any infected individual may transmit the disease to a susceptible one according to the following relation:

where γ denotes the per capita recovery rate. I 0 is the initial condition for the infected population.

• Once an individual has been moved from Infected to Recovered, it is assumed that it is not possible to be infected again. This condition is described by:

where R 0 is the initial condition for the recovered population.

It is important to emphasize that the population size (N) along time t is defined as

In practice, the model parameters must be determined to represent a particular epidemic. For this purpose, it is necessary to formulate and to solve an inverse problem. In the section that describes the methodologies adopted in this work, more details on the formulation and solution of this problem is presented.

Mathematically, an OCP can be formulated as follows (Bryson and Ho, 1975; Feehery and Barton, 1996; Lobato, 2004) . Initially, let

where z is the vector of state variables, and u is the vector of control variables. Ψ and L are the first and second terms of the performance index, respectively. The minimization problem is given by arg min

with consistent initial conditions given by

According to the optimal control theory (Bryson and Ho, 1975; Feehery and Barton, 1996) , the solution of the OCP, whose problem is defined by Eqs. (5) and (6), is satisfied by the co-state equations and the stationary condition given, respectively, byλ

where H is the Hamiltonian function defined by

This system of equations is known as the Euler-Lagrange equations (optimality conditions), which are characterized as Boundary Value Problems (BVPs). Thus, to solve this model, an appropriated methodology must be used, as for example, Shooting Method or Collocation Method (Bryson and Ho, 1975) . As mentioned by Bryson and Ho (1975) and Feehery and Barton (1996) , the main difficulties associated with OCPs are the following: the existence of end-point conditions (or region constraints) implies multipliers and associated complementary conditions that significantly increase the complexity of solving the BVP using a indirect method; the existence of constraints involving the state variables and the application of slack variables method may introduce differential algebraic equations of higher index; the Lagrange multipliers may be very sensitive to the initial conditions.

Differential Evolution is a powerful optimization technique to solve mono-objective optimization problems, proposed by Storn and Price (1997) . This evolutionary strategy differs from other population-based algorithms in the schemes considered to generate a new candidate to solution of the optimization problem (Storn and Price, 1997; Price et al., 2005; Lobato and Steffen, 2011) . The population evolution proposed by DE follows three fundamental steps: mutation, crossover and selection. The optimization process starts by creating a vector containing NP individuals, called initial population, which are randomly distributed over the entire search space. During G max generations, each of the individuals that constitute the current population are subject to the procedures performed by the genetic operators of the algorithm.

In the first step, the mutation operator creates a trial vector by adding the balanced difference between two individuals in a third member of the population, by v (G+1)

The parameter F represents the amplification factor, which controls the contribution added by the vector difference, such that F ∈ [0, 2]. In turn, Storn and Price (1997) proposed various mutation schemes for the generation of trial vectors (candidate solutions) by combining the vectors that are randomly chosen from the current population, such as:

The second step of the algorithm is the crossover procedure. This genetic operator creates new candidates by combining the attributes of the individuals of the original population with those resulting in the mutation step. The vector u (G+1)

jk , such as k = 1, . . . , d, where d denotes the dimension of the problem and randb (k) ∈ [0, 1] is a random real number with uniform distribution. The choice of the attributes of a given individual is defined by the crossover coefficient, represented by CR, such that CR ∈ [0, 1] is a constant parameter defined by the user. In turn, rnbr ( j) ∈ [1, d] is a randomly chosen index.

After the generation of the trial vector by the steps of mutation and crossover, the evolution of the best individuals is defined according to a greedy strategy, during the selection step. Price et al. (2005) have defined some simple rules for choosing the key parameters of DE for general applications. Typically, one might choose NP in the range from 5 to 10 times the dimension (d) of the problem. In the case of F, it is suggested taking a value ranging between 0.4 and 1.0. Initially, F = 0.5 may be a good choice. In the case of premature convergence, F and NP might be increased.

The multi-objective optimization problem (MOP) is an extension of the mono-objective optimization problem. Due to the conflict between the objectives, there is no single point capable of optimizing all functions simultaneously. Instead, the best solutions that can be obtained are called optimal Pareto solutions, which form the Pareto curve (Deb, 2001) . The notion of optimality in a MOP is different from the one regarding optimization problems with a single objective. The most common idea about multi-objective optimization found in the literature was originally proposed by Edgeworth (1881) and further generalized by Pareto (1896) . One solution is said to be dominant over another, if it is not worse in any of the objectives, and if it is strictly better in at least one of the objectives. As an optimal Pareto solution dominates any other feasible point in the search space, all of these solutions are considered better than any other. Therefore, multi-objective optimization consists of finding a set of points that represents the best balance in relation to minimizing all objectives simultaneously, that is, a collection of solutions that relates the objectives, which are in conflict with each other, in most cases.

Let f (x) = ( f 1 (x) , . . . , f m (x)) T be the objective vector such that f k : P → IR, for k = 1, . . . , m, where x ∈ P is called decision vector and its entries are called decision variables. Mathematically, a MOP is defined as (Deb, 2001; Lobato, 2008) :

Due to the favorable outcome of DE in solving mono-objective optimization problems, for different fields of science and engineering, Lobato and Steffen (2011) proposed the Multi-Objective Differential Evolution (MODE) algorithm to solve multi-objective optimization problems. Basically, this evolutionary strategy differs from other algorithms by the incorporation of two operators to the original DE algorithm, the mechanisms of rank ordering (Deb, 2001; Zitzler and Thiele, 1999) , and exploration of the neighborhood for potential solution candidates (Hu et al., 2005) . A brief description of the algorithm is presented next.

At first, an initial population of size NP is randomly generated, and all objectives are evaluated. All dominated solutions are removed from the population by using the operator Fast Non-Dominated Sorting (Deb, 2001) . This procedure is repeated until each candidate vector becomes a member of a front. Three parents generated by using DE algorithm are selected at random in the population. Then, an offspring is generated from these parents (this process continues until NP children are generated). Starting from population P 1 of size 2NP, neighbors are generated to each one of the individuals of the population. These neighbors are classified according to the dominance criterion, and only the non-dominated neighbors (P 2 ) are put together with P 1 , in order to form P 3 . The population P 3 is then classified according to the dominance criterion. If the number of individuals of the population P 3 is larger than a predefined number, the population is truncated according to the Crowding Distance (Deb, 2001) criterion. This metric describes the density of candidate solutions surrounding an arbitrary vector. A complete description of MODE is presented by Lobato and Steffen (2011) .

As mentioned earlier, the first objective of this work is to determine the parameters of the SIR model adopted to predict the evolution of COVID-19 epidemic considering experimental data from China. In this case, it is necessary to formulate and to solve an inverse problem. It arises from the requirement of determining parameters of theoretical models in such a way that it can be employed to simulate the behavior of the system for different operating conditions. Basically, the estimation procedure consists of obtaining the model parameters by the minimization of the difference between calculated and experimental values.

In this work, it is assumed that, since the outbreak persists for a relatively short period of time, the rate of births and deaths of the population is insignificant. Thus, we take µ = 0, since there are probably few births/deaths in the corresponding period. We are interested in the determination of the following parameters of the SIR model: β, γ and I 0 . Let and I sim i are the experimental and simulated infected population, respectively, and M represents the total number of experimental data available. In this case, the SIR model must be simulated considering the parameters calculated by DE, in order to obtain the number of infected people estimated by the model and, consequently, the value of the objective function (F ). As the number of measured data, M, is usually much larger than the number of parameters to be estimated, the inverse problem is formulated as a finite dimensional optimization problem in which we aim at minimizing F .

In order to formulate both OCPs, the parameters estimated considering the proposed inverse problem are used. As proposed by Neilan and Lenhart (2010) and Biswas et al. (2014) , a new variable W, which denotes the number of vaccines used, is introduced in order to determine the optimal control strategy for vaccine administration. For this purpose, the total amount of vaccines available during the whole period of time is proportional to uS . Physically, u represents the portion of susceptible individuals being vaccinated per unit of time (Biswas et al., 2014) . It is important to mention that u acts as the control variable of such system. If u is equal to zero there is no vaccination, and u equals to one indicates that all susceptible population is vaccinated. A schematic diagram of the disease transmission among the individuals for the SIR model with vaccination is shown in Fig. 2 .

Mathematically, the SIR model considering the presence of control is written as: 

where W 0 is the initial condition for the total amount of vaccines. It is important to emphasize that the population size (N) after the inclusion of this new variable W along the time t is defined as N(t) = S (t) + I(t) + R(t) + W(t).

The first formulation aims to determine the optimal vaccine administration (u) to minimize the infected population, represented by Ω 1 . Thus, let

The OCP is defined as arg min

u Ω 1 (16) subject to Eqs. (11) -(14) and u min ≤ u ≤ u max , where t 0 and t f represents the initial and final time, respectively, and u min and u max are the lower and upper bounds for the control variable, respectively. The second formulation considers two objectives, i.e., the determination of the optimal vaccine administration, in order to minimize the number of infected individuals and, at the same time, to minimize the number of vaccines needed. The total number of vaccines can be determined by

whereas the number of infected people is given by Eq. (15) . Thus, the multi-objective optimization problem is formulated as arg min u (Ω 1 , Ω 2 ) (18) subject to Eqs. (11) -(14) and u min ≤ u ≤ u max . In both problems, the control variable u must be discretized. In this context, the approach proposed consists of transforming the original OCP into a nonlinear optimization problem. For this purpose, let the time interval 0, t f be discretized using N elem time nodes, with each node denoted by t i , where i = 0, . . . , N elem − 1, such that t 0 ≤ t i ≤ t f . For each of the N elem − 1 subintervals of time, given by [t i , t i+1 ], the control variable is considered constant by parts, that is,

In order to obtain an optimal control strategy for vaccine administration, that can be used in medical practice, we consider the bang-bang control which consisting of a binary feedback control that turns either "on" (in our case, when u = u max = 1) or "off" (when u = u min = 0) at different time points, determined by the system feedback. In this case, as the control strategy u is constant by parts, the proposed optimal control problem has N elem − 2 unknown parameters, since the control variable at the start and end times are known. The resulting nonlinear optimization problems are solved by using the DE, in the case of the mono-objective problem, given by Eq. (16), and MODE, for the multi-objective problem defined by Eq. (18).

In order to apply the proposed methodology to solve the inverse problem described previously, the following steps are established:

• Objective function: minimize the functional F , given by Eq. (10);

• Design space: 0.1 ≤ β ≤ 0.6, 0.04 ≤ γ ≤ 0.6 and 10 −8 ≤ I 0 ≤ 0.5 (all defined after preliminary executions);

• DE parameters: population size (25), number of generations (100), perturbation rate (0.8), crossover rate (0.8) and strategy rand/1 (as presented in Section 4.1). The evolutionary process is halted when a prescribed number of generations is reached (in this case, 100). 20 independent runs of the algorithm were made, with different seeds for the generation of the initial population;

• To evaluate the SIR model during the optimization process, the Runge-Kutta-Fehelberg method was used;

• Initial conditions: S (0) = 1 − I 0 , I(0) = I 0 , and R(0) = 0. In this case, I 0 is chosen as the first reported data in relation to the number of infected individuals in the time series;

• The data used in the formulation of the inverse problem refer to the population of China, from January 22 to April 2, 2020, taken from Johns Hopkins Resource Center (2020 (accessed April 03, 2020). Table 1 presents the results (best and standard deviation) obtained using DE. It is possible to observe that DE was able to obtain good estimates for the unknown parameters and, consequently, for the objective function, as can be verified, by visual inspection of Fig. 3 . These results were obtained, as mentioned earlier, from 20 runs. Thus, the values of the standard deviation demonstrate that the algorithm converges, practically, to the same optimum in all executions (best). Physically, the probability of transmission by contact in the Chinese population is superior to 35 % (β equal to 0.3566). In addition, γ equal to 0.0858 implies a moderate per capita recovery rate. One must consider that, since many cases may not be reported, for different reasons, as for example an asymptomatic infected person, the value of I 0 may vary, as well as the behavior of the model over time. It is important to emphasize that when choosing I 0 as a design variable, the initial condition for the susceptible population (S 0 ) is automatically defined, that is, S 0 = 1 − I 0 , since there is not, at the beginning of an epidemic, a considerable number of recovered individuals and, thus, R 0 = 0 is a reasonable choice. In this case, the available data refer to the number of infected individuals and these represent only the portion of individuals in the population that have actually been diagnosed. This is due, among other facts, to the lack of tests to diagnose the disease of all individuals who present symptoms. Thus, as the number of susceptible individuals at the beginning of the epidemic is dependent on the value of I 0 , in this work it is considered that the total size of the population, typically defined as N = S + I + R, is actually a portion of the total population, since the number of infected individuals available is also a fraction of those who have actually been diagnosed. In this case, the results presented below represent only the fraction of the infected population that was diagnosed and, consequently, the fraction of individuals susceptible to contracting the disease. Qualitatively, the results presented are proportional to the number of individuals in the population who were diagnosed with the disease.

In order to evaluate the sensitivity of the solutions obtained, in terms of the objective function, the best solution (β = 0.3566, γ = 0.0858, and I 0 = 0.0038) was analyzed considering a perturbation rate given by δ. For this purpose, the range [(1 − δ) θ k , (1 + δ) θ k ] was adopted, for k ⊂ {1, 2, 3}, where θ = (β, γ, I 0 ). Thus, in each analysis, one design variable is perturbed and the value of F in relation to this noise is computed. Figure 4 presents the sensitivity analysis for each estimated parameter, in terms of the objective function, considering δ equal to 0.25 and 100 equally spaced points in the interval of interest. In these figures, it is possible to observe that the variation in the result of each parameter, as expected, in a worst value for the F . In addition, that the design variable more sensible to δ parameter is the β parameter, since a wide range of values for the F were obtained. 

We consider two distinct analysis in this section, in order to evaluate the proposed methodology considered to solve the mono-objective optimization problem: i) solution of the proposed mono-objective optimal control problem and ii) evaluation on the influence of the maximum amount of vaccine, by defining an inequality constraint. For this purpose, the following steps are established:

• Objective function: minimize the functional Ω 1 , given by Eq. (16);

• The previously calculated parameters (β, γ and I 0 ) are employed in the simulation of the SIR model;

• Design space: 0 ≤ t i ≤ t f , for i = 1, . . . , t N elem −1 , and N elem = 10. It is important to mention that this value was chosen after preliminary runs, i.e., increasing this value do not produce better results in terms of the objective function;

• DE parameters: population size (25), number of generations (100), perturbation rate (0.8), crossover rate (0.8) and strategy rand/1 (as presented in Section 4.1). The evolutionary process is halted when a prescribed number of generations is reached (in this case, 100). 20 independent runs of the algorithm were made, with different seeds for the generation of the initial population;

• To evaluate the SIR model during the optimization process, the Runge-Kutta-Fehelberg method was used;

• Initial conditions: S (0) = 1 − I 0 , I(0) = I 0 , and R(0) = 0. As in the previous case, I 0 is chosen as the first reported data in relation to the number of infected individuals in the time series; Table 2 presents the best solution obtained by using DE and considering ten control elements, in terms of the number of individuals. The objective function obtained (about 8945.4278 individuals) is less than the case in which no control is considered (about 1594607.2234 individuals), i.e., the number of infected individuals is lower when a control strategy is considered (see Figs. 5(a) and 5(c)). If the number of infected individuals is reduced, due to control action, the number of susceptible individuals rapidly decreases until its minimum value (1.4382 × 10 −3 individuals) and, consequently, the number of recovered individuals rapidly increase until its maximum value (767.5187 individuals), as observed in Figs. 5(b) and 5(d), respectively. In terms of the action regarding the control variable, the effectiveness is readily verified in the beginning of the vaccine administration. Further the administration is conducted in specific intervals of time, which preserves the health of the population, as observed in Fig. 5(e) . The evolution of the number of vaccinated individuals is presented in Fig. 5(f) . In this case, due to control action, the vaccinated population increase rapidly until the value is saturated (141835.1405 individuals). In summary, all obtained profiles are coherent from the physical point of view. Finally, it is important to mention that the standard deviation for each result is, approximately, equal to 10 −3 , which demonstrates the robustness of DE to solve the proposed mono-objective optimal control problem. In this model, the evaluation of the number of vaccinated individuals is associated with an inequality constraint. This relation bounds the quantity of individuals that can be vaccinated due to the limitation related to the production of vaccines. For this purpose, two control elements are incorporated to the model: if W(t 1 ) ≤ W lim , then u = 1. Otherwise, u = 0 (t 1 is the instant of time that W(t 1 ) = W lim , and W lim is the upper bound for the number of vaccinated individuals). Table 3 presents the results obtained considering different quantities for the parameter W lim . As expected, the insertion of this constraint implies in limiting the maximum number of vaccinated individuals and, consequently, a lower number of individuals are vaccinated. The increase of the parameter W lim implies in the reduction of the objective function value, in number of infected and recovered individuals and, consequently, an increase in the number of susceptible individuals. These analysis can be observed in Fig. 6 . 

As presented previously, a multi-objective optimal control problem was proposed in order to minimize the number of infected individuals (Ω 1 ) and to minimize the quantity of vaccine administered (Ω 2 ). To evaluate the proposed methodology considered to solve this multi-objective optimization problem, the following steps are established:

• Objective functions: minimize both Ω 1 and Ω 2 together, which are defined by Eqs. (15) and (17), respectively;

• The previously calculated parameters (β, γ and I 0 ) are employed in the simulation of the SIR model;

• Design space: 0 ≤ t i ≤ t f , for i = 1, . . . , t N elem −1 , and N elem = 10. • MODE parameters: population size (50), number of generations (100), perturbation rate (0.8), crossover rate (0.8), number of pseud-curves (10), reduction rate (0.9), and strategy rand/1 (as presented in Section 4.1). The stopping criterion adopted is the same as in the previous cases.

• To evaluate the SIR model during the optimization process, the Runge-Kutta-Fehelberg method was used;

• Initial conditions: S (0) = 1 − I 0 , I(0) = I 0 , R(0) = 0, and W(0) = 0. Table 4 . It must be stressed that the Pareto curve presents the non-dominated solutions, as described in Section 4.2. The point A represents the best solution in terms of the minimization of the number of infected individuals, with Ω 1 = 8963.7775, that is, the number of infected individuals at t f assume its lowest value, which is equal to I(t f ) = 769.0921, but considering a larger amount of vaccine administered (Ω 2 = 6.9358). On the other hand, the point B represents the best solution in terms of the quantity of vaccine administered, with Ω 2 = 1.2940, i.e, the minimization of such value when t = t f . However, for this point, the number of infected individuals is high (Ω 1 = 56644.0350). The point C is a compromise solution, which is a good solution in terms of both objectives simultaneously, with intermediary values for both objectives -Ω 1 = 13298.2440 and Ω 2 =2.3034. Table 4 . In Figure 7 (e) it is possible to observe the activation of the control variable when vaccine is introduced. Besides, in both results obtained, the action of such treatment is readily verified in the population during a larger interval of time in the beginning of the vaccine administration. In Figures 7(b) , 7(c), 7(d) and 7(f) the susceptible, infectious, recovered and number of vaccines profiles are presented, respectively, for each point described in Table 4 . In these figures we can visualize the importance of the control strategy used. For example, the points A and C are good choices in terms of the minimization of infected individuals, although the point A has a highest value in terms of the objective Ω 2 . On the other hand, point B is satisfactory in terms of minimizing the amount of vaccines administered, but, from a clinical point of view, it is not a good choice, as the number of infected individuals is not minimized.

In this contribution was proposed and solved an inverse problem to simulate the dynamic behavior of novel coronavirus disease (COVID-19) considering real data from China. The parameters of the compartmental SIR (Susceptible, Infectious and Recovered) model were determined by using Differential Evolution (DE). Considering the parameters obtained with the solution of the proposed inverse problem, two optimal control problems were proposed. The first consists of minimizing the quantity of infected individuals. In this case, an inequality that represents the quantity of vaccines available was analyzed. The second optimal control problem considers the minimizing together the quantity of infected individuals and the prescribed vaccine concentration during the treatment. This problem was solved using Multi-Objective Differential Evolution (MODE). In general, the solution of the proposed multi-objective optimal control problem provides information from which an optimal strategy for vaccine administration can be defined.

The use of mathematical models associated with optimization tools may contribute to decision making in situations of this type. It is important to emphasized that the quality of the results is dependent on experimental data considered. In this context, one may cite the following limitations regarding the SIR model: i) poor quality of reported official data and ii) the simplifications of the model, as for example terms as birth rate, differential vaccination rate, weather changes and its effect on the epidemiology. Finally, it is worth mentioning that the problem formulated in this work is not normally considered in the specialized literature (only the minimization of the infected individuals is normally proposed). In this context, the formulation of the multi-objective optimization problem and its solution by using MODE represents the main contribution of this work. 

Modeling and Analysis of an SEIR Epidemic Model with a Limited Resource for Treatment

Numerical Modeling and Theoretical Analysis of a Nonlinear Advection-reaction Epidemic System

Scheduling of Measles Vaccination in Low-income Countries: Projections of a Dynamic Model

Application of the ARIMA Model on the COVID-2019 Epidemic Dataset

Advances in Simultaneous Strategies for Dynamic Process Optimization

A SEIR Model for Control of Infectious Diseases with Constraints

An Introduction to Compartmental Modeling for the Budding Infectious Disease Modeler

Parameter Estimation Based Synchronization for an Epidemic Model with Application to Tuberculosis in Cameroon

Applied Optimal Control: Optimization, Estimation and Control

A Familial Cluster of Pneumonia Associated with the 2019 Novel Coronavirus Indicating Person-to-person Transmission: a Study of a Family Cluster

Forecasting the Spread of Mosquito-Borne Disease using Publicly Accessible Data: A Case Study in Chikungunya

Multi-objective Optimization Using Evolutionary Algorithms

Deep Learning to Detect Alzheimer's Disease from Neuroimaging: a Systematic Literature Review

Dynamic Simulation and Optimization with Inequality Path Constraints

Predicting Unobserved Exposures from Seasonal Epidemic Data

Severe Acute Respiratory Syndrome-related Coronavirus: the Species and Its Viruses-A Statement of the Coronavirus Study Group

The mathematics of infectious diseases

A New Multi-objective Evolutionary Algorithm: Neighbourhood Exploring Evolution Strategy

Clinical Features of Patients Infected with 2019 Novel Coronavirus in Wuhan

Mapping 2019-nCoV

Modeling Infectious Diseases in Humans and Animals

Global Stability of an SEIR Epidemic Model Where Empirical Distribution of Incubation Period is Approximated by Coxian Distribution

Dynamic Analysis of an SEIR Model with Distinct Incidence for Exposed and Infectives

Global Stability for the SEIR Model in Epidemiology

Early Transmission Dynamics in Wuhan, China, of Novel CoronavirusâȂŞInfected Pneumonia

A Conceptual Model for the Coronavirus Disease 2019 (COVID-19) Outbreak in Wuhan, China with Individual Reaction and Governmental Action

Hybrid Approach for Dynamic Optimization Problems

Multi-objective Optimization for Engineering System Design

Determination of an Optimal Control Strategy for Drug Administration in Tumor Treatment Using Multi-objective Optimization Differential Evolution

A New Multi-objective Optimization Algorithm Based on Differential Evolution and Neighborhood Exploring Evolution Strategy

Developing Covid-19 Vaccines at Pandemic Speed

Mathematical Analysis of a Sex-structured HIV/AIDS Model with a Discrete Time Delay

An introduction to optimal control with an application in disease modeling

Cours d'Économie Politique. F. Rouge

Modelling the Effect of Changes in Vaccine Effectiveness and Transmission Contact Rates on Pertussis Epidemiology

The Effect of Control Strategies to Reduce Social Mixing on Outcomes of the COVID-19 Epidemic in Wuhan, China: a Modelling Study

Differential Evolution: A Practical Approach to Global Optimization

Pattern of Early Human-to-human Transmission of Wuhan 2019 Novel Coronavirus (2019-nCoV)

Why is it Difficult to Accurately Predict the COVID-19 Epidemic? Infectious Disease Modelling

Inference and Forecast of the Current West African Ebola Outbreak in Guinea

Stability of SEIR Model of Infectious Diseases with Human Immunity

Differential Evolution-A Simple and Efficient Heuristic for Global Optimization over Continuous Spaces

Deterministic SEIRs Epidemic Model for Modeling Vital Dynamics, Vaccinations, and Temporary Immunity

Clinical Characteristics of 138 Hospitalized Patients With 2019 Novel CoronavirusâȂŞInfected Pneumonia in Wuhan, China

Radiotherapy Workflow and Protection Procedures During the Coronavirus Disease 2019 (COVID-19) Outbreak: Experience of the Hubei Cancer Hospital in Wuhan

The SIR Model and the Foundations of Public Health

Susceptible Exposed Infected Recovery (SEIR) Model with Immigration: Equilibria Points and its Application

The Epidemiology, Diagnosis and Treatment of COVID-19

Multiobjective Evolutionary Algorithms: a Comparative Case Study and the Strength Pareto Approach

This study was financed in part by the Coordenação de Aperfeiçoamente de Pessoal de Nível Superior-Brasil (CAPES)-Finance Code 001, Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (FAPERJ), and Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq).