Determination of an optimal control strategy for vaccine administration in COVID-19 pandemic treatment.

BACKGROUND AND
OBJECTIVE: For decades, mathematical models have been used to predict the behavior of physical and biological systems, as well as to define strategies aiming at 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 the novel coronavirus disease (COVID-19) is considered an important theme due to the quantity of infected people worldwide. In this work, the objective is to determine an optimal control strategy for vaccine administration in COVID-19 pandemic treatment considering real data from China. Two optimal control problems (mono- and multi-objective) to determine a 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.
METHODS: An inverse problem is formulated and solved in order to determine the parameters of the compartmental Susceptible-Infectious-Removed model. The solutions for both optimal control problems proposed are obtained by using Differential Evolution and Multi-objective Optimization Differential Evolution algorithms.
RESULTS: A comparative analysis on the influence related to the inclusion of a control strategy in the population subject to the epidemic is carried out, in terms of the compartmental model and its control parameters. The results regarding the proposed optimal control problems provide information from which an optimal strategy for vaccine administration can be defined.
CONCLUSIONS: The solution of the optimal control problem can provide information about the effect of vaccination of a population in the face of an epidemic, as well as essential elements for decision making in the economic and governmental spheres.

Introduction

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, and promotion of the improvement of various health-economic aspects [1]. 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 [2], [3]. One can cite several studies using models for measles vaccination [4], [5], HIV/AIDS [6], tuberculosis [7], dengue [8], pertussis epidemiology [9], 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 November 2019, that has rapidly spread around in many other countries worldwide [10], [11]. The common symptoms are severe respiratory illness, fever, cough, and myalgia or fatigue, especially at the onset of illness [12]. The transmission may happen person-to-person, through direct contact or droplets [13], [14], [15].

Since the COVID-19 outbreak in Wuhan city in November of 2019, various computational model-based predictions have been proposed and studied. Lin et al. [16] 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 period of time. Benvenuto et al. [17] proposed the use of an Auto Regressive Integrated Moving Average (ARIMA) model to predict the spread, prevalence and incidence of COVID-19. Roda et al. [18] used a Susceptible-Infectious-Removed (SIR) model to predict the COVID-19 epidemic in Wuhan after the lockdown and quarantine. In this study, the authors demonstrate that non-identifiability in model calibrations using the confirmed-case data is the main reason for wide variations in the results. Prem et al. [19] 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 aiming 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 [20] studied a SEIR model with nonlinear incidence rates in epidemiology, in terms of global stability of endemic equilibrium. Al-Sheikh [1] evaluated a SEIR epidemic model with limited resources for treating infected people. For this purpose, the existence and stability of disease-free and endemic equilibrium were investigated. Li and Cui [21] 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. [22] 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. [5] proposed a SEIR model with immigration and determined the system equilibrium conditions. Kim et al. [23] developed a Coxian-distributed SEIR model considering an empirical incubation period, and a stability analysis was also performed.

In order to reduce the dissemination of COVID-19 worldwide, various procedures have been adopted. As mentioned by Zhai et al. [24] and Wei et al. [25], 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. To date, according to the U.S. Food and Drug Administration, there are no specific therapies to COVID-19. To our best knowledge, however, treatments including antiviral agents, corticosteroids, antibodies, convalescent plasma transfusion and radiotherapy are being studied [26].

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 to COVID-19, but there is a huge effort to develop a vaccine in a record time, which justifies the present study [27]. Mathematically, the determination of optimal protocols for vaccine administration characterizes an optimal control problem. This particular optimization problem consists in the determination of control variable profiles that minimize (or maximize) a given performance index [28], [29]. In order to solve this problem, several numerical methods have been proposed [28], [30], [31], [32]. These methods are classified according to three broad categories: direct optimization methods, Pontryagin’s Maximum Principle based methods and Hamilton-Jacob-Bellman based methods. The direct approach is the most traditional strategy considered to solve an optimal control problem, 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 [30].

From an epidemiological point of view, Neilan and Lenhart [33] 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. [34] studied different mathematical formulations for an optimal control problem considering a Susceptible-Exposed-Infectious-Removed model. For this purpose, these authors analyzed 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 imposing 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 work, 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 characterize the proposed mathematical model (based on the compartmental SIR model), an inverse problem is formulated and solved using the Differential Evolution algorithm [35], [36]. 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. In order to solve each problem, Differential Evolution and Multi-objective Optimization Differential Evolution algorithms [37] are employed, respectively.

This work is organized as follows. In Section 2, the fundamental details of the compartmental model are shown, as well as the general formulation of optimal control problems, and the metaheuristics used to solve the proposed problems. This section also presents the inverse problem used to retrieve the parameters of the epidemiological model, in addition to the optimal control problems proposed (mono- and multi-objective). The results obtained by solving the inverse problem and the optimal control problems are presented in Section 3. In turn, discussions on the results are presented in Section 4. Finally, conclusions are drawn in Section 5.

Methods

Mathematical modeling in epidemiology

In the specialized literature, various compartmental models used to represent the evolution of an epidemic can be found [9], [38], [39], [40], [41]. The study of these models is very important to understand the epidemic spreading mechanisms and, consequently, to investigate the transmission dynamics in population [38]. As mentioned by Keeling and Rohani [42], these compartmental models can be divided into two groups: i) population-based models and; ii) agent-based 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 integro-differential 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 yet 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 I—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 classical compartmental models can be represented as [42], [43]:

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

Susceptible-Infected-Removed (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-Removed (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 [42] and Hethcote [43], 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 and; 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 COVID-19 epidemic in China. The choice of this model is due to the study conducted by Roda et al. [18]. 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 .

Fig. 1

Compartments in the SIR model [42].

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 . These groups (S, Iand R) can be scaled by population size N using the normalized variables S n, I n and R n, defined as S n ≡ S/N, I n ≡ I/N and R n ≡ R/N. Thus the population is normalized in such a way that and we have the new system: where S n0, I n0 and R n0 are the initial condition for the susceptible, infected and recovered populations, respectively. These normalized variables denote the fractions of the number of individuals in classes S, I and R in relation to population N. In practice, the model parameters must be determined to represent a particular epidemic. For this purpose, the formulation of an inverse problem is presented in Section 2.2, which is solved using the Differential Evolution algorithm, briefly described in Section 2.4. In turn, Eqs. (4)–(6) are solved using the Runge-Kutta-Fehlberg technique, presented in Section 2.3.

Formulation of the inverse problem

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 reported data from China. In this case, it is necessary to formulate and 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 reported 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 by natural cases or other reasons of the population is insignificant. Thus, we take since there are probably few births/deaths in the corresponding period. We are interested in the determination of the following parameters of the normalized SIR model: β, γ and I n0. It is important to mention that I n0 is used to define the initial condition of all dependent variables of the model. Thus, consider that the number of infected people, obtained from the solution of Eqs. (4)–(6), is stored in a vector represented by I, and let the function

Mathematically, the inverse problem is formulated as subject to Eqs. (4)–(6), where and I are the reported and simulated infected population in normalized form for the i-th unit of time, respectively, is the highest reported value for the infected normalized population, and M represents the total number of reported data available. In this case, the normalized SIR model must be simulated considering the parameters calculated by Differential Evolution, in order to obtain the number of infected people estimated by the model and, consequently, the value of the objective function (). 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 [44].

Runge-Kutta-Fehlberg technique for ordinary differential equations

The Runge-Kutta-Fehlberg technique [45] is an extension of the classical Runge-Kutta methods for numerical solution of ordinary differential equations. The technique aims to increase the accuracy of Runge-Kutta methods by adopting a strategy of step size adaptation. Essentially, at each step, the method calculates the difference (absolute error) between the fourth and fifth order approximations of the solution. If the resulting difference is below an arbitrary threshold, the calculated value is accepted. Otherwise, the step size is reduced and new approximations are calculated. On the other hand, in case the two approximations are in very close agreement, the step size is increased. For the sake of brevity, the method is not detailed here, but only its main attributes. More details are presented by Mathews and Fink [46].

Let Υ be an open set of  ×  . Consider an initial value problem given by with f: Υ ⊂  ×   →  and called the initial condition, where (t 0, y 0) ∈ Υ. Taking the step size h, each step of the Runge-Kutta-Fehlberg method requires calculating the values

When y(t) is known, the fourth order approximation to y(t) at is given byIn turn, the fifth order approximation is given by Eqs.  (9) and (10) are iteratively compared, in order to decide on adjusting the step size h.

Differential evolution

Differential Evolution (DE) is a stochastic population-based metaheuristic, proposed by Storn and Price [35], to solve single-objective optimization problems over continuous spaces. Essentially, the evolutionary strategy follows three fundamental steps: mutation, crossover and selection. The initial population containing NP individuals is randomly created, covering the entire search space. The population in a given generation G is composed of d-dimensional individuals denoted by for . During G maxgenerations, the three genetic operators are applied sequentially, so that one hopes the population evolves towards the optimizer of the problem.

In the first step, mutant vectors are created by adding the balanced difference between two individuals to a third individual, by means of . Individuals are mutually different and selected at random, and Frepresents the scale factor, such that F ∈ [0, 2]. In the crossover procedure, the second step, new candidate solutions are created by combining the attributes of the original population with the mutant vectors. Thus, trial vectors are created by if randb(k) ≤ CR or . Otherwise, where and randb(k) ∈ [0, 1] is an uniformly distributed random number. The crossover probability CR ∈ [0, 1] is a predefined constant parameter. In turn, rnbr(j) ∈ [1, d] is a randomly chosen index. After generating the trial vectors, the best individuals are selected according to a greedy strategy, during the third step. Further details on the algorithm can be found in Price et al. [36].

Optimal control problems

Consider a system whose fundamental property is to operate in several states. Let x(t) ∈  be the vector of state variables of the system at time t ∈ [0, t f], with t f > 0. In turn, the vector of control variables at time t is denoted by u(t) ∈  . Given the initial state the evolution of the system is described by the state equation where f:  ×   ×  →  is a continuously differentiable vector function, and the notation is usually used to represent dx(t)/dt.

If the initial state and the control trajectory—the set of values of u(t) over the interval [0, t f]—are known, the result of the integration of Eq.  (11) provides the state trajectory x(t) over the same interval. Given the continuously differentiable functions F:  ×   ×  → and ϕ:  ×  → . In an optimal control problem, the objective is to calculate u(t) so that the state and control trajectories minimize the objective function while subject to boundary conditions .

In general, control variables may be constrained according to u(t) ∈ Ω(t), where Ω(t) ⊂  is the set of feasible values at time t. Instead, optimal control problems may also be constrained by mixed inequality constraints given by g(x(t), u(t), t) ≤ 0, and constraints involving only state variables, such as h(x(t), t) ≤ 0, where g:  ×   ×  →  and h:  ×  →  are continuously differentiable functions for each t ∈ [0, t f] [47].

According to Bryson Jr and Ho [28], in order to find the control trajectory that produces a stationary value of the objective given by Eq.  (12), first consider the Hamiltonian where H:  ×   ×   ×  → . The solution of the general optimal control problem is obtained by solving both the adjoint equation and the state equation given by Eq.  (11), where u(t) is obtained from . The boundary conditions are given by .

This system of equations is known as the Euler-Lagrange equations, which defines necessary conditions for optimality. In order to solve this model, an appropriated methodology must be used, as for example, the Shooting Method or the Collocation Method. The main difficulties associated with optimal control problems are the following: (i) multipliers and associated complementary conditions can significantly increase the complexity of solving the problem using an indirect method due to the existence of end-point conditions (or region constraints); (ii) the Lagrange multipliers may be very sensitive to the initial conditions and; (iii) constraints involving the state variables and the application of slack variables method may introduce differential algebraic equations of higher index [28]. In the next section we propose an optimal control problem for vaccine administration, based on the compartmental model presented in Section 2.1.

Optimal control problem for vaccine administration

In order to formulate both optimal control problems, the parameters estimated considering the proposed inverse problem are used. As proposed by Neilan and Lenhart [33] and Biswas et al. [34], a new variable W, which denotes the number of vaccinated individuals, 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 [34]. 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 vaccination is taking place and all susceptible population will be vaccinated as time goes towards to infinity. A schematic diagram of the disease transmission among the individuals for the normalized SIR model with vaccination is shown in Fig. 2 .

Fig. 2

Compartments in the normalized SIR model with vaccination.

Mathematically, the normalized SIR model considering the presence of control (referred to here as SIRW) is written as: where W n0 is the initial condition for the total amount of vaccinated individuals in normalized form. It is important to emphasize that after the inclusion of this new variable W n along the time t, the normalized form in relation to population is defined as .

The first formulation aims to determine the optimal vaccination policy to minimize a quantity related to the normalized infected population, represented by Ω1. Thus, let

The optimal control problem is defined as subject to Eqs.  (13)–(16) and u min ≤ u ≤ u max, where t 0 and t f represents the initial and the 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 vaccination policy, in order to minimize Ω1 and, at the same time, to minimize some metric related to the number of vaccines needed, represented by Ω2. For this purpose, we opted by the following objective function

It should become clear that other quantities could be employed for this second objective, such as ∫dW n. Thus, the multi-objective optimization problem is formulated as subject to Eqs.  (13)–(16) and u min ≤ u ≤ u max. In both problems, the control variable u must be discretized. In this context, the approach proposed consists on transforming the original optimal control problem 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, where such that t 0 ≤ t ≤ t f. For each of the subintervals of time, given by the control variable is considered constant by parts, that is, for where u min ≤ u ≤ u max.

In order to obtain an optimal control strategy for vaccination policy, that can be used in medical practice, we consider the bang-bang control which consists of a binary feedback control that turns either “on” (in our case, when ) or “off” (when ) at different time points, determined by the system feedback. In this case, as the control strategy uis constant by parts, the proposed optimal control problem has unknown parameters, since the control variable at the start and end times are known. The mono-objective optimization problem, given by Eq.  (18), is solved using Differential Evolution, presented in Section 2.4. In turn, the basic concepts regarding multi-objective optimization are presented in Section 2.7, and the problem defined by Eq.  (20) is solved using Multi-objective Optimization Differential Evolution, described in Section 2.8.

Multi-objective optimization problem

Let denote a hyperparallelepiped of all . The vector is called decision vector, its entries are the decision variables, and it is bounded by for every . The domain is known as the decision variables space. In turn, assume that with m ≥ 2, is called the multi-objective function, where and for are the objectives. Throughout this work, the vector x is referred to as a point or a solution interchangeably, since it is within the set of elements in the decision variables space that eventually represents a solution for the optimization problem.

In addition, consider the functions and . The set of equality and inequality constraints that delimits the subspace to be searched for the optimal solution are expressed as g(x) ≤ 0and where g(x) ≤ 0, for and for . The set of solutions that satisfy both equality and inequality constraints is known as the feasible setor search space, and denoted by . A point that belongs to Ω is termed as a feasible solution. In turn, the objective space, represented by is the image of the multi-objective function f.

The multi-objective problem is defined by

Here, we focus on minimization problems. For maximization problems, the concepts presented here must be adjusted accordingly. In general, there is no single point capable of minimizing all functions simultaneously in a multi-objective problem, due to the conflict between the objective functions. It means that none of the objectives can be improved without deterioration to at least one of the other objectives. This is attained by the concepts of dominance of solutions and Pareto optimality, which are presented below.

The concept of domination defines a way to compare potential solutions in optimization problems with multiple objectives. A solution x 1 is said to dominate another solution x 2, denoted by x 1 ≺ x 2, if both of the conditions below are satisfied:In other words, a 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. Solution x 1 is said not to dominate solution x 2 in case any of these conditions is violated.

Solution x 1 is no worse than x 2 in all objectives, that is, f(x 1) ≤ f(x 2), for all ;

Solution x 1 is strictly better than x 2 in at least one objective, that is, f(x 1) < f(x 2) for some .

Multi-objective optimization consists in finding a set of feasible solutions, the Pareto optimizers, which in turn are denoted by x* ∈ Ω, that represents the best balance regarding the minimization of all objectives simultaneously. An objective vector z* ∈ V is said to be a Pareto optimum if there is no other objective vector z ∈ V such that for every and for some index . Thus, the so-called Pareto set, denoted by Θand also referred to as Pareto front or Pareto curve (the latter holds for ), usually has an infinite number of Pareto optimal values. The objective vector z* is a Pareto optimum if . Moreover, .

Multi-objective optimization differential evolution

Due to the favorable outcome of DE in solving mono-objective optimization problems in different fields of science and engineering, Lobato and Steffen Jr [37] proposed the Multi-objective Optimization Differential Evolution (MODE) algorithm to solve multi-objective optimization problems. The fundamental structure of MODE is composed of the computational procedures briefly described below. Initially, a population of N individuals is randomly generated in the decision variables space. Dominated candidates are removed from the population, and the remaining points are ranked in several Pareto fronts according to the FastNonDominatedSort procedure proposed by Deb et al. [48]. Next, three candidate solutions are selected at random. From this group of individuals, a new candidate solutions is generated and incorporated in the population. The process of selecting individuals and generating new candidates is carried out repeatedly until N new points are produced, causing the population to reach 2N individuals.

In the next step, the algorithm searches for dominant points in the neighborhood of each point in the current population using the neighborhood exploring evolution strategy, proposed by Hu et al. [49]. The purpose of this procedure is to create new candidate solutions around the 2N current points, by the action of small perturbations in the population. In this way, newly generated offsprings replace their dominated parents, accelerating the convergence. Then the population is truncated according to the CrowdingDistance criterion defined by Deb et al. [48]. This step is intended to select the best individuals, in order to reduce the population size back to N for the next generation. This iterative procedure goes until a maximum number of generations is reached. At the end of the generations, it is expected that MODE has converged on an approximate Pareto set, with adequate diversity of solutions. A complete description of MODE is presented by Lobato and Steffen Jr [37].

Results

This section presents the parameters retrieved by solving the proposed inverse problem using data related to the number of infected individuals in China, from January 22 to April 2, 2020, taken from Ref. [50]. These parameters are used to solve optimal control problems, in order to assess the impact of vaccination of the population during the epidemic.

In all simulations presented here, single-objective optimization problems are solved using Differential Evolution with 25 individuals in the population, perturbation rate and crossover rate equal to 0.8. In the case of multi-objective optimization problems, the solutions are obtained using Multi-objective Optimization Differential Evolution with population size equal to 50, perturbation rate and crossover rate equal to 0.8, reduction rate equal to 0.9, and 10 pseudo-curves. For both metaheuristics, the iterative procedure is halted when a prescribed number of generations is reached (in this case, 100). In both cases, the control parameters were independently adjusted and tests with populations containing more individuals showed no significant gain in terms of results.

Inverse problem

In order to obtain the parameters of the SIR model that best fit the time series of infected individuals, the inverse problem given by Eq.  (8) is minimized. The design space, defined after preliminary executions [16], [19], is bounded by 0.1 ≤ β ≤ 0.6, 0.04 ≤ γ ≤ 0.6, and . The initial conditions for the normalized system of ordinary differential equations are and .

Table 1 presents the results (best and standard deviation) calculated from 20 independent runs of the optimization algorithm. It is possible to observe that Differential Evolution 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 for the dimensional variables S, I and R. Thus, the values of the standard deviation demonstrate that the algorithm converges, nearly, to the same optimum in all executions (best).

Table 1

Results obtained for the proposed inverse problem considering Differential Evolution.

	β(day^−1)	γ(day^−1)	In0	F
Best	0.3566	0.0858	0.0038	0.1649
Standard Deviation	1.2545×10^−5	1.6291×10^−5	1.43238×10^−6	1.2260×10^−7

Fig. 3

Simulated and reported profiles considering the estimated parameters.

Considering the data of Table 1, we can obtain an estimate for the value of the basic reproduction number calculated by . We can observe that this value is coherent with that recently reported by Sanche et al. [51]. These authors found for COVID-19 in Wuhan, which is quite superior from previous estimates [52]. As pointed out by Wearing et al. [53], we may expect differences in the estimation of from i) initial epidemic growth rate and ii) trajectory matching (this last one used in our approach). Keeling and Grenfell [54] also discuss how different approaches–using deterministic or stochastic modelling–affect the calculation of .

In order to evaluate the sensitivity of the solutions obtained, in terms of the objective function, the best solution was analyzed considering a perturbation rate given by δ. For this purpose, the range was adopted, for , where . Thus, in each analysis, one design variable is perturbed and the value of in relation to this noise is computed. Fig. 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.

Fig. 4

Sensitivity analysis of estimated parameters.

Mono-objective optimal control problem

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, consider the objective function Ω1, given by Eq.  (18). The previously calculated parameters (β, γ and I n0) are employed in the simulation of the normalized SIRW model, with initial conditions given by and . For this problem the design space is bounded by 0 ≤ t ≤ t f, for and . It is important to mention that this value was chosen after preliminary runs and increasing N elemdo not produce better results in terms of the objective function. Table 2 presents the best solution obtained in terms of the number of individuals in dimensional form.

Table 2

Results obtained for the proposed mono-objective optimization problem in dimensional form (days).

Ω1 (Number of Individuals  ×  Days)	S(tf)	I(tf)	R(tf)	W(tf)
8945.4278	1.4382×10^−3	2.1201	767.5187	141835.1405

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 . Otherwise, (t 1 is the instant of time that 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. Results related to this approach are shown in Fig. 6 (in dimensional form).

Table 3

Results obtained for the proposed mono-objective optimization problem considering different quantities for the parameter Wlim (days).

Wlim	t1 (Days)	Ω1 (Number of Individuals  ×  Days)	S(tf)	I(tf)	R(tf)
50000	1.4389	907790.2114	9674.7066	6659.1695	77888.4004
100000	2.1985	76039.6424	35197.4174	1472.4917	6524.2013

Fig. 6

Influence of the maximum amount of vaccine in the objective function, susceptible-infectious-removed populations profiles, control variable strategy and number of vaccinated individuals’ profiles (in dimensional form).

Multi-objective optimal control problem

As presented previously, a multi-objective optimal control problem was proposed in order to simultaneously minimize Ω1 and Ω2. Thus, the optimization problem represented by Eq.  (20) is solved in order to obtain an optimal Pareto set, considering the previously calculated parameters β, γ and I n0 for the normalized SIRW model, where the initial conditions and the design space are the same as in Section 3.2.

Fig. 7 (a) presents the Pareto curve and three points (A, B and C) belonging to this curve, as shown in Table 4 . It must be stressed that the Pareto curve presents the non-dominated solutions, as described in Section 2.8. In Fig. 7(b)–7(f) are presented the susceptible-infectious-removed populations profiles (in dimensional form), control variable strategy and number of vaccinated individuals’ profiles (in dimensional form) considering the points presented in Table 4.

Fig. 7

Pareto curve, susceptible-infectious-removed populations profiles, control variable strategy and number of vaccinated individuals’ profiles (in dimensional form).

Table 4

Some points belonging to the Pareto curve obtained by proposed multi-objective optimization problem (in dimensional form) (Days).

Point	Ω1 (Number of Individuals  ×  Days)	Ω2	S(tf)	I(tf)	R(tf)	W(tf)
A	8963.7775	6.9358	135.2256	2.1653	769.0921	141698.2905
B	56644.0350	1.2940	33749.2312	907.8714	4860.0582	103087.6259
C	13298.2440	2.3034	13697.0991	20.4643	1140.9894	127746.2276

Discussion

According to the results obtained by solving the inverse problem, 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. Regarding the sensitivity analysis of estimated parameters, shown in Fig. 4, it is possible to observe that, as expected, the variation leads to a deteriorated value, in relation to the optimum value obtained by the optimizer. In addition, the most sensitive parameter, that is, that causes the greatest variation in the objective value, is β. Clearly, such uncertainty can affect the prediction of the model behavior.

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, 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 it represents 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 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 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 the mono-objective optimal control problem, the objective function obtained (about 8945.4) is less than the case in which no control is considered (about 1594607.2), 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 () and, consequently, the number of recovered individuals rapidly increase until its maximum value (767.5 individuals), as observed in Fig. 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). 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 which demonstrates the robustness of Differential Evolution to solve the proposed mono-objective optimal control problem.

Fig. 5

Objective function, susceptible-infectious-removed populations profiles in dimensional form, control variable strategy and number of vaccinated individuals’ profiles in dimensional form.

As expected, considering a constraint related to the maximum amount of vaccine 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, as can be seen in Fig. 6.

Now, consider the multi-objective problem. The point A, in Fig. 7(a), represents the best solution in terms of the minimization of the first objective, with but considering a larger amount of vaccine administered . On the other hand, the point B represents the best solution in terms of the second objective, with i.e, the minimization of such value when . However, for this point, the attainment of the first objective is poor (). The point C is a compromise solution, which is a good solution in terms of both objectives simultaneously, with intermediary values for both objectives, and . It must be stressed that the choice of the points A, B and C is absolutely arbitrary and merely illustrative of the compromise between the objectives Ω1and Ω2.

Fig. 7 (e) shows 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 Fig. 7(b), 7(c), 7(d) and 7(f) the susceptible, infectious, recovered and number of vaccinated individuals 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 first objective, 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 second objective, but from a clinical point of view, it is not a good choice, as the metric related to the number of infected individuals is not minimized.

Conclusions

In this contribution it is 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. Considering the parameters obtained with the solution of the proposed inverse problem, two optimal control problems were proposed. The first consists on minimizing a quantity related to the number of infected individuals. In this case, an inequality that represents the quantity of vaccines available was analyzed. The second optimal control problem considers minimizing the same metric used in the mono-objective problem together a second quantity related to the number of vaccines used during the treatment. This problem was solved using Multi-objective Optimization Differential Evolution. 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 emphasize that the quality of the results is dependent on the quality of the reported data considered. In this context, one may cite the following limitations regarding the SIR model: ) poor quality of reported official data and; ) 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 Multi-objective Optimization Differential Evolution represents the main contribution of this work.

Ethical approval

The authors state that ethical approval is not required for this type of study in the country where the study was conducted.

Declaration of Competing Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be constructed as a potential conflict of interest. The authors have no affiliation with any organization with a direct or indirect financial interest in the subject matter discussed in the manuscript.