Optimal strategy for a dose-escalation vaccination against COVID-19 in refugee camps

1.   Introduction

The world is in the midst of the COVID-19 (coronavirus disease 2019) pandemic caused by the new virus severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2). As of November 30, 2020, over 61 million cases of COVID-19 have been confirmed, accounting for 1.5 million deaths globally ^[1]. More seriously, the COVID-19 epicentre has transferred to developing countries, such as Russia, Brazil and India, with rapidly increasing infection levels. Especially, the challenge of second wave would increase in the winter. Thus, more refugee camps have been exposed to the risk of COVID-19, with even more dire consequences than observed in general populations ^[2].

Due to the unavailability of both specific medications to treat the disease and extensive screening, COVID-19 has a high transmission rate, especially in refugee camps with poor access to water and limited sanitation. A safe and effective vaccine against SARS-CoV-2 would obviate the current need for social distancing, which is impractical for those who reside in settlements with high density. Fortunately, there are more than 100 COVID-19 vaccine candidates under development ^[3]. In particular, a phase 2 trial on the immunogenicity and safety of COVID-19 vaccines in healthy adults has been closed, and a vaccine is believed to be available soon ^[4,5,6,7,8]. These findings provide a foundation for identifying optimal vaccine administration strategies. However, in the short term, the vaccine supply could not cover the worldwide need. To facilitate equitable access and distribution of these vaccines ^[9], the criteria for COVID-19 vaccine prioritization are a focus of attention. Hence, it is worth studying optimal vaccine allocation and administration under limited resources.

For decades, the role of mathematical models in evaluating the spread and control of infectious diseases cannot be overemphasized. These models can also proffer solutions for transmission control, treatment and vaccination by employing infection assumptions and mathematical approaches. Accurate outbreak prediction and setting of efficient control strategies significantly depend on high-quality mathematical models. As an epidemic transmission prediction method, Kermack and McKendrick proposed the classic compartment model of infectious disease dynamics in 1927, which divides the population into compartments ^[10]. Thereafter, differential equations and kinetic systems methods have become the mainstream approaches for quantitative description of epidemic transmission. For COVID-19, the Susceptible-Infectious-Recovered (SIR) model and its variants ^[11,12,13] have been used to characterize transmission in the early stage. Soon afterwards, the Susceptible-Exposed-Infectious-Recovered (SEIR) model with the latent compartment has become a consensus model among studies of SARS-CoV and has provided the basis for model settings for COVID-19 epidemiological dynamics studies ^[14]. To date, studies based on the Susceptible-Exposed-Asymptomatic-Infectious-Recovered (SEAIR) model have become mainstream ^[2,15,16]. Clearly, the model selection is updated as medical researchers gain deeper insights into the COVID-19 transmission mechanism. Numerous models implying detection status, quarantine situations, and asymptomatic/presymptomatic infections have boomed, with the aim of understanding the mechanisms underlying COVID-19 spread and to project transmission dynamics of the outbreak ^{[17,18,19,20,21]}. Thus, it is possible for us to study COVID-19 from a rigorous mathematical perspective by using the theory of differential equations. Specifically, according to the stability of equilibrium points, we can analyse whether COVID-19 will become endemic ^[22]. And based on sensitivity analysis, efficient control methods can be employed ^[23], which provides promising primary directions for governments to fight against the epidemic.

However, it should be noted that many parameters are involved in compartment models, and the number of parameters increases as more complex transmission mechanisms are considered. To make the model in accordance with historical data, overfitting would occur, resulting in poor prediction ability. Hence, many researchers turn to artificial intelligence (AI) for help. For example, machine learning models, which are represented by deep learning models due to their extraordinary approximation capability, were established to predict the outbreak size in China and throughout the world in ^[24,25], respectively. Other methods that work directly with empirical data, such as the Bayesian method ^[26], Prophet forecasting procedure ^[27], and autoregressive integrated moving average method ^[28], have also been developed. The above modelling methods can achieve acceptable accuracy, and neither compartment division nor parameter estimation procedures are required, which is an advantage over traditional compartment model-based methods. We cannot ignore the significant contributions of AI techniques against COVID-19 (not only limited to the modelling processes discussed above but also including detection & diagnosis ^[29,30,31], infodemiology & infoveillance ^[32,33,34], and biomedicine & pharmacotherapy ^[35,36,37]). However, clear transmission mechanisms are not reflected in such models. In our opinion, it would be more appealing to develop AI-assisted modelling methods that consider fusion with physical laws.

At the global level, to fight against the spread of COVID-19, control measures have been implemented by governments. Currently, in the absence of an effective and safe vaccine, control strategies primarily rely on home quarantine, public traffic bans, lockdowns and social distancing ^[21,38]. Mathematical models have been shown to provide factual guides in proffering mitigation measures for infectious diseases. More importantly, they can also be used to optimize containment and intervention measures, balancing cost-effectiveness and infection control. Especially under limited resources, optimal control techniques play a significant role in infectious disease intervention and have been widely applied in the Ebola, SARS, MERS and HPV epidemics ^[39,40]. Recently, an increasing number of studies have explored optimization leading to the best outcomes at the population outbreak level, taking into consideration the economic constraints concerning COVID-19. Among these studies, ^[41] solved an optimal control problem of minimizing deaths and control measure implementation costs using the forward-backward sweep method while waiting for the vaccine. A similar numerical method was also applied in ^[42], considering hospitalization rates and environmental spraying. In ^[22], an optimal control problem was formulated introducing governmental intervention measures and taking into account economic constraints. High-quality control problems containing various factors of contact rate and isolation rate using a multi-objective genetic algorithm have been proposed ^[43]. Along with the above studies in the context of open-loop optimal control, closed-loop control techniques have also been implemented in recent work. For example, ^[44] studied the optimal screening and testing level using the model predictive control (MPC) technique. Taking model uncertainty into account, ^[45] further developed a robust MPC-based feedback policy using interval arithmetic that adapts social distancing measures cautiously and safely. Advanced control techniques can greatly contribute to the fight against COVID-19.

Optimal control theory offers a reliable approach to optimize a given objective in a nonlinear, dynamic system. Treatments, including antiviral agents, antibodies and radiotherapy, are still ambiguous, and greater expectations have been placed on the efficiency of vaccines against SARS-CoV-2. Vaccines are being developing in record time; mathematically, the optimal protocol for vaccination characterizes an optimal control problem. Meanwhile, studies on the optimal control of vaccination are still very few. In ^[46], researchers computed the threshold of a hypothetical imperfect vaccine that could lead to elimination of COVID-19 in the United States. In ^[47], a multi-objective optimal control problem was solved considering minimization of vaccine concentrations based on the SIR compartmental model using single/multi-objective differential algorithms. The global COVID-19 pandemic is now severe, and this condition is considered to remain for the foreseeable future. Medical resources (including medicine, equipment and health workers) are still in seriously short supply, especially for regions such as refugee camps. However, in most studies on optimal control of COVID-19, only extremely simple constraints are imposed on the control variable (i.e., the control variable is simply restrained within the interval [0, 1]), but practical constraints on limited resources cannot be mirrored faithfully. Hence, optimal vaccine administration under limited medical resources is worth further study.

In this paper, we study the optimal strategy of a dose-escalation vaccination approach for COVID-19 considering limited medical resources. The SEAIRD compartment model is taken as the basic model to characterize the transmission dynamics. Two vaccination strategies, i.e., the single-dose and double-dose strategies, are considered by extending the basic model. Taking the vaccination rates in susceptible compartments as control variables, the issue is solved under the framework of nonlinear constrained optimal control problems. Along with the box constraint imposed on the control variable, two extra kinds of practical constraints on medical resources are considered. On the one hand, the total vaccine supply is restrained. On the other hand, considering the shortage of medical care personnel, the maximum number of individuals vaccinated daily is restrained. Taking the Kutupalong-Balukhali refugee camp in Bangladesh as the simulation background, the optimal control performances of the two vaccination strategies are simulated and then compared based on various indices. The results suggest that both strategies can effectively reduce COVID-19 dissemination, while the double-dose strategy would achieve better control. Further studies reveal that increasing the vaccine supply and providing more medical care personnel can greatly contribute to epidemic control. The main contribution can be summarized into the following four items:

(1) The optimal vaccine administration strategy for COVID-19 is solved under the framework of nonlinear optimal control, meeting the balance between infection number and vaccine consumption.

(2) Necessary constraints on medical resources are considered in the formulated optimal control problem, making the computed vaccine administration strategy more practical.

(3) The optimal control performance under two scenarios, i.e., single-dose and double-dose vaccination administration, are considered and evaluated.

(4) We reveal that better COVID-19 control can be achieved by increasing the vaccine supply and the number of medical care personnel.

The remainder of the paper is organized as follows: Section 2 describes the formulation of the mathematical model representing the evolution of the COVID-19 pandemic in the setting of refugee camps. In Section 3, we present the optimal vaccination control problem. Section 4 analyses the characterization of the optimal control strategy. In Section 5, numerical simulations are carried out, and two vaccination strategies are compared and discussed from various aspects. Finally, conclusions and potential future research directions are summarized in Section 6.

2.   Formulation of the epidemic model

The mathematical model describes the high COVID-19 transmission rate in the setting of refugee camps is formulated in this section. Our model is adapted from the model proposed in ^[2] with the key parameters and initial values referred to in ^[2,16]. In Section 2.1, we introduce the classic SEIR model. In Section 2.2, a modified SEAIRD model is established to fit the COVID-19 outbreak in refugee camps. Then in Sections 2.3 and 2.4, the established SEAIRD model is extended to describe one-dose and double-dose vaccine administration, respectively.

2.1.   The basic SEIR model

In the context of compartment models, in general, the total population $N$ can be divided into several typical compartments, e.g., $S$ (susceptible, not yet infected), $E$ (exposed, being infected, in the incubation period, no symptoms), $I$ (infected, infectious) and $R$ (recovered, no longer infectious). Considering the spreading pattern, the compartmental models can be represented as follows: Susceptible-infected (SI), susceptible-infected-removed (SIR), susceptible-infectious-susceptible (SIS), susceptible-exposed-infectious-removed (SEIR), etc.

The classical SEIR model divides the population into four compartments: susceptible ($S$), exposed ($E$), infected ($I$) and recovered ($R$). Unlike the SIR model, the SEIR model further considers the latent period in which a virus carrier without immediate illness is capable of infecting other susceptible individuals. The presence of the latent period results in a longer pandemic transmission cycle.

The total population $N$ is fixed in the classical SEIR model, which assumes that the virus spreads within a closed system. The individuals in the four compartments are homogeneously mixed, and the individuals move between compartments through a certain probability of interaction. The infected $I$ would turn the susceptible $S$ into the exposed $E$ with the SARS-CoV-2 transmission rate $\beta$; $\alpha$ represents the portion of the exposed $E$ that would become the infected $I$; and the infected individuals $I$ become cured with a certain probability $\gamma$ and turn into the recovered $R$. Meanwhile, $N = S+E+I+R$ is satisfied. The set of differential equations is as follow:

2.2.   An improved SEAIRD model

The dynamic characteristic of COVID-19 transmission fitted to refugee camps is based on the SEIR model and considers asymptomatic individuals with reference to the potential COVID-19 outbreak in the Kutupalong-Balukhali refugee camp in Bangladesh ^[2] and models describing the Italian COVID-19 epidemic ^[16]. The total population is divided into six compartments, i.e., susceptible (could fall ill through contact with infected individuals) $S$, exposed (incubating, not infectious, no symptoms) $E$, asymptomatic infected (infected, infectious, no symptoms/presymptomatic) $A$, actively infected (infected, infectious, with symptoms) $I$, recovered (healing) $R$ and dead (extinct) $D$. Considering the extremely limited number of hospital beds in refugee camps, quarantine is not taken into account. Hence, the transmission dynamics can be established as

where ${N}_{0}\left(t\right) = S\left(t\right)+E\left(t\right)+A\left(t\right)+I\left(t\right)+R\left(t\right)$. Individuals in $S$ proceed into the exposed compartment $E$ at a rate $\beta$ through contact with infected individuals in $A$ or $I$, and then, they reside for an incubation period. Considering the lack of healthcare capacity and great survival pressure, the transmission rate is similar for asymptomatic infected individuals and symptomatic individuals. Exposed individuals become infected at a rate $\alpha$, where a fraction $\delta$ experiences symptomatic infection and the remainder experience asymptomatic infection. The recovery times for actively infected and asymptomatic infected patients are given by $1/{\gamma }_{I}$ and $1/{\gamma }_{A}$, respectively. The disease-induced mortality rate in the $A$ compartment is denoted as $\sigma$.

Such an SEAIRD compartmental model has been widely used to describe the evolution of COVID-19. Some preliminary results have been conducted by previous researches. Here, we present some important relevant results directly, and interesting readers can refer to ^[48]. By using the next generation method, the basic production number of model (2.2) can be derived as follow:

2.3.   Introducing single-dose vaccine administration into the SEAIRD model

In this paper, we assume that vaccines will be available to fight against SARS-CoV-2, and the effectiveness of a single-dose vaccine ${\nu }_{1}$ is estimated according to trial phase 1/2 of the ChAdOx1 nCoV-19 vaccine ^[6].

where ${N}_{1}\left(t\right) = {S}_{0}\left(t\right)+{S}_{1}\left(t\right)+E\left(t\right)+A\left(t\right)+I\left(t\right)+R\left(t\right)$.

The corresponding compartment diagram is illustrated in Figure 1, where ${S}_{0}$ denotes susceptible individuals without vaccination and ${S}_{1}$ represents susceptible individuals with one-dose vaccination. Let the individuals in ${S}_{0}$ be vaccinated at a rate of ${w}_{1}$. The vaccinated individuals become immune at a ratio of ${\nu }_{1}$, proceeding into the $R$ compartment. If neutralizing antibody responses against SARS-CoV-2 are not detected ^[5], the risk of the vaccinated individuals in ${S}_{1}$ exposed to the infection is assumed as well as the general susceptible individuals in ${S}_{0}$, and the identical transmission rate $\beta$ is used.

2.4.   Introducing double-dose vaccine administration into the SEAIRD model

A trial report of the ChAdOx1 nCoV-19 vaccine in the UK showed that anti-spike IgG responses rose after a booster dose ^[6]. Hence, we characterize how a double-dose vaccine administration policy would influence the dynamics of the COVID-19 epidemic as follow:

where ${N}_{2}\left(t\right) = {S}_{0}\left(t\right)+{S}_{1}\left(t\right)+{S}_{2}\left(t\right)+E\left(t\right)+A\left(t\right)+I\left(t\right)+R\left(t\right)$.

The corresponding compartment diagram is illustrated in Figure 2, where ${S}_{0}$ denotes susceptible individuals without vaccination, ${S}_{1}$ and ${S}_{2}$ represent susceptible individuals administered the first and the second vaccine dose, respectively. ${S}_{0}$ would be first vaccinated at a rate of ${u}_{1}$, and ${u}_{2}$ for ${S}_{1}$ with the second dose. ${\nu }_{1}$ and ${\nu }_{2}$ represent the effectiveness of the first and second vaccine dose, respectively.

3.   Establishing optimal control problems

3.1.   Single-dose vaccine administration strategy

We introduce the variable $V$ to record the number of vaccinated individuals in compartment ${S}_{0}$ as follow:

It is assumed that the maximum number of daily vaccinated individuals is proportional to the number of medical care personnel. The following three types of constraints are considered when constructing the optimal control problem:

(1) Considering the limited medical care personnel, the maximum number of individuals vaccinated daily is constrained as follow:

(2) The total vaccine supply is considered to be limited, which can be expressed as the following constraint:

(3) The vaccination rate ${w}_{1}$ naturally exists in the range of $\left[\mathrm{0, 1}\right]$, i.e.,

With the above settings, we can formulate the optimal single-dose vaccine administration problem as the following nonlinear constrained optimal control problem, which is denoted as Problem 1:

3.2.   Double-dose vaccine administration strategy

We introduce the following variable to record the number of vaccinated individuals in compartments ${S}_{0}$ and ${S}_{1}$ :

Following the ideas in Section 3.1, the following constraints are taken into account in this scenario:

Thus, the optimal double-dose vaccine administration problem, which is denoted as Problem 2, can be formulated as follow:

4.   Characterization of optimal control

In this section, we characterize the structure of the optimal control problems proposed in Section 3 by the parametric variational principle. Since Problem 2 is more complicated, only its characterization is derived in detail. Before we derive the following results, it is necessary to guarantee existence and uniqueness of the solution to Problem 2, and related proofs are provided in Appendix A.

First, the inequality constraints in Problem 2 are transformed into equality constraints by non-negative parametric variables as follows:

where $\boldsymbol{\omega } = {\left[{\omega }_{1}, {\omega }_{2}, {\omega }_{3}, {\omega }_{4}, {\omega }_{5}, {\omega }_{6}\right]}^{T}$ is the parametric variable list. Thus, according to ^[46], the Hamiltonian of Problem 2 is obtained as follow:

where $\boldsymbol{\lambda } = {\left[{\lambda }_{{S}_{0}}, {\lambda }_{{S}_{1}}, {\lambda }_{{S}_{2}}, {\lambda }_{E}, {\lambda }_{A}, {\lambda }_{I}, {\lambda }_{R}, {\lambda }_{V}\right]}^{T}$ is the adjoint variable (or the costate variable in the context of optimal control) that introduces system equations into the Hamiltonian and $\boldsymbol{\mu } = {\left[{\mu }_{1}, {\mu }_{2}, {\mu }_{3}, {\mu }_{4}, {\mu }_{5}, {\mu }_{6}\right]}^{T}$ is the non-negative parametric variable list that introduces equality constraints into the Hamiltonian. Applying the parametric variational principle ^[49] to Problem 2, for the optimal solutions, the following first-order necessary conditions must be satisfied:

with the following transversality conditions:

where Eq (4.3) is equivalent to the system equation; Eq (4.4) is called the adjoint equation, which can be further expressed as

where $K = \frac{\beta \left(A+I\right)}{{N}_{2}^{2}}$, and the time variable $t$ is omitted for simplicity.

The optimal control and corresponding state trajectory are denoted as ${\boldsymbol{u}}^{*} = {\left[{u}_{1}^{*}, {u}_{2}^{*}\right]}^{T}$ and ${\boldsymbol{x}}^{*} = {\left[{S}_{0}^{*}, {S}_{1}^{*}, {S}_{2}^{*}, {E}^{*}, {A}^{*}, {I}^{*}, {R}^{*}, {V}^{*}\right]}^{T}$, respectively. By solving Eq (4.5), the optimal control is determined as follows:

In most epidemic control problems where only box constraints are imposed on control variables (i.e., Eqs (3.9) and (3.10)), when arriving at equations similar to Eqs (4.16) and (4.17), one can further characterize optimal control as min-max formulations where no parametric variable is involved by using the standard argument for control bounds ^[50]. Though an extra mixed state-control constraint (i.e., Eq (3.7)) and pure-state constraint (i.e., Eq (3.8)) are also involved in Problem 1, one can follow the analysis in ^[51] to obtain similar characterization, i.e.,

Then, together with the system equations and the adjoint equations, the forward-backward sweep method ^[50] that is commonly used in epidemic control can be used to solve the problem.

However, two control variables are coupled in the mixed state-control constraint in Problem 2. Hence, it is not possible to further simplify Eqs (4.16) and (4.17), which suggests that the forward-backward sweep method cannot guarantee numerical stability herein. Actually, the forward-backward sweep method is only applicable to optimal control problems with fixed initial state boundary conditions and free terminal state boundary conditions. The characterization of optimal control in a simple formulation must be provided to implement the numerical process. However, such limitations do not exist in advanced computational optimal control techniques, such as direct collocation methods ^[52,53] and indirect symplectic pseudospectral methods ^{[49,51,54,55]}.

5.   Numerical simulations and discussion

In this section, we evaluate the effectiveness of the two vaccine administration strategies. The optimal control problem was solved using the ICLOCS ^[52] software package. In detail, the Hermite method with 1000 nodes was used to implement the discretization, and the IPOPT software package ^[56] was used to solve the resultant nonlinear programming problem. By using ICLOCS, under the framework of direct transcription method, the optimal control problem is transcribed into a nonlinear programming (NLP) in the following fashion:

where $\boldsymbol{X}$, as the unknown quantities to be determined, consists of state and control variable at the discretization points.

5.1.   Parameter settings

The parameters used for numerical simulations mainly refer to the literature ^{[2,6,16,57,58]} and are listed in Table 1. The effectiveness of vaccines is based on the phase 1/2 trial report of the ChAdOx1 nCoV-19 vaccine ^[6]. The initial total population was set as $600000$. In both scenarios, we set ${S}_{0} = 599988$, ${E}_{0} = 10$ and ${I}_{0} = {A}_{0} = 1$, while the initial populations in other compartments are zero. The final time is set as ${t}_{f} = 90$.

5.2.   Comparison between two vaccination strategies

We explore the control effects of two vaccination strategies in this sub-section. Observing the cost functional in both Problems 1 and 2, they are both constituted of two parts, i.e., one term related to infected individuals and another term related to control intensity. In the following simulations, we set ${b}_{1} = {b}_{2} = 0.001$ in both two vaccination strategies. As in the single-dose scenario, we set $d = 1$; and ${c}_{1} = {c}_{2} = 1$ are selected in the double-dose scenario. The optimal state trajectory of two scenarios are reported in Figure 3. And the profiles of corresponding optimal control are given in Figures 4 and 5. The cost functions obtained in the single-dose and the double-dose scenarios are 979.26 and 581.99, respectively.

From Figure 3, it is seen that both two vaccination strategies can effectively control the spread the COVID-19. To give a better comparison between these two strategies, we summarize their control effects in Table 2. Additionally, the profiles of accumulated vaccinated individuals and accumulated death are reported in Figure 6. One can read that the double-dose vaccination strategy was found to further lower the peaks of $A$ and $I$ compared with those in the single-dose vaccination strategy. A smaller number of deaths occurred with the double-dose vaccination strategy than with the other strategy, which is approximately only 10% of the uncontrolled case. The population in $A$ at ${t}_{f}$ in the single-dose vaccination scenario was much higher than that in the uncontrolled scenario. Moreover, we note that available vaccines are not used up in the single-dose scenario. It is mainly because that, when compared with the double-dose scenario, more individuals in the susceptible compartments become infected and then proceed into other compartments. The above comparisons suggest that the double-dose vaccination strategy can achieve better COVID-19 control than the single-dose strategy. However, we also noted that the single-dose vaccination strategy postpones the times corresponding to the peaks of $A$ and $I$ more efficiently. More susceptible individuals could be vaccinated under the single-dose vaccination strategy. In contrast, fewer individuals would gain immunity limited to vaccine supply under the double-dose vaccination strategy. Subsequently, individuals with a single-dose vaccine might also be infected, and thus, the peak of infection and death would be higher but come later than with the double-dose vaccination scenario.

One may get confused that at the final time, population of the $R$ compartment in the single-dose vaccination is slightly lower than that in the uncontrolled scenario. It can be drawn from Figures 4 and 5 that vaccination stops at approximately the 65th day. It suggests that from that day on, the simulation is almost a free-spread process. When we look back at the system dynamics in Eqs (2.4) and (2.5), without vaccination, only infected individuals (individuals in the $A$ or $I$ compartments) can proceed into the $R$ compartment by recovery. From Figure 3(d), (e), it is seen that there're more infected individuals in the uncontrolled scenarios when compared with the controlled scenario. Hence, the population in the $R$ compartment at a much higher speed.

5.3.   Effects of daily vaccination capabilities and total vaccine supply

Figures 4(b) and 5(b) show that in both vaccination strategies, the number of daily vaccinations arrives at a maximum level in the initial period. This partially suggests that further increasing the number of medical care personnel, which results in a higher daily vaccination capability $\Omega$, would lead to better control performance. Moreover, noting that the supplied vaccines are exhausted in the double-dose vaccination strategy, we may infer that further increasing the total vaccine supply ${V}_{max}$ would reach a more satisfactory control performance. To validate the above two conjectures, in Figure 7, we test the effect of $\Omega$ and ${V}_{max}$ on several indices in the double-dose scenario.

The results remind us that increasing the vaccine supply and equipping more medical care personnel are critical in the fight against COVID-19. If the single-dose and double-dose vaccination simulations are combined, a more feasible strategy could be proposed. In the early stage, more susceptible individuals are supposed to be vaccinated with the first dose to postpone the infection peak and thus buy time for vaccine production. A double-dose vaccination strategy should be implemented at a later period. In this way, the objectives of both the number and peak of infections could be balanced.

6.   Conclusions and future research directions

6.1.   Conclusions

Two vaccine administration strategies, i.e., single-dose and double-dose strategies, for COVID-19 in refugee camps considering limited medical resources were studied in this paper based on extended SEAIRD models. Numerical simulations show that both strategies can control the spread and efficiently reduce the total infections by integrating the optimal control technique. Interestingly, the double-dose strategy performs better on infection demonstration than the signal-dose strategy under the same vaccine constraints. This result implies that standard vaccine administration with a double dose for a portion of the population should be satisfied, rather than more populations being vaccinated with one dose first. However, if additional control measures are implemented or medical resources are replenished, postponing the peak of infection would be more significant. For the challenging settings of refugee camps, vaccines alone are not enough; epidemic control also depends on the maximum level of daily vaccination, which demonstrates that additional health workers are as important as vaccines themselves.

As the simulations show, the vaccines cannot entirely stop the spread of COVID-19, and discovery of effective drugs for treatment is still necessary and expected. Moreover, the trial results for different vaccine types and different trial phases are inconsistent, and SARS-CoV-2 has been mutating ^[59]. Thus, the availability of a vaccine could decline. We make optimistic presumptions in this work, while the challenges in the real world are more severe.

6.2.   Future research directions

(1) Recent research has revealed that the successful immune rate severely depends on age ^[7]. Hence, it seems interesting to consider the age-structure in future work. One simple idea is to further divide the ${S}_{0}$, ${S}_{1}$ and ${S}_{2}$ compartments into more sub-compartments according to the age-structure. Another interesting idea is to describe the epidemic dynamics with partial differential equations (PDEs), where the coefficients are functions of the age variable. Correspondingly, efficient numerical techniques to solve optimal control problems for PDEs are required therein.

(2) Vaccination strategies are formulated over a relatively long time span. Thus, it seems more practical to set the parameters in the model to be time-dependent ^[51]. For example, the successful immune rate will generally decrease over time considering drug resistance. In addition, the incidence rate would vary seasonally.

(3) Minimal numerical analysis is given in this paper. In the following works, numerical analysis such as the stability, parameter sensitivity, second derivative of Lyapunov and strength number ^[20] will be supplemented.

(4) Only one control measure, i.e., vaccination, is considered in this paper. In future work, compound control strategies should be studied to attain more satisfactory control performance.

Appendix A: Existence and uniqueness of the solution to Problem 2

Problem 2 can be rewritten into the following fashion:

subject to system dynamics and initial conditions

and state and control constraints

Compared with Eq (3.11), the total vaccine supply constraints are transformed into an equivalent form on the terminal value of state $V$.

The boundedness of system (A.1) for the finite time interval is proved first. We note that the populations in compartments ${S}_{0}$, ${S}_{1}$, ${S}_{2}$, $E$, $A$ and $I$ decreases proportionally to their current quantities, it guarantees that all these variables remain non-negative as their initial values are non-negative. As for state variables $V$ and $R$, their initial values and derivatives are non-negative, thus their non-negativity can be guaranteed. To establish the upper bounds for state variables, we consider the total population size ${N}_{2}$. The change in ${N}_{2}$ satisfies that ${\dot{N}}_{2}\left(t\right) = -\alpha I\left(t\right)$ and is ${N}_{2}$ bounded by its initial condition ${N}_{2}\left(0\right)$. Since ${S}_{0}$, ${S}_{1}$, ${S}_{2}$, $E$, $A$, $I$ and $R$ are all non-negative, the upper bound of ${N}_{2}$ is also their upper bound. As for the upper bound for state variable $V$, it allows from the boundedness of control variables ${u}_{1}$ and ${u}_{2}$ as well as state variables ${S}_{1}$ and ${S}_{2}$.

We give a generalized formulation of the problem (A.1)–(A.3) as follow:

subject to system dynamics

initial boundary conditions

terminal boundary conditions

and constraints

Assume that the functions $F:{\mathbb{R}}^{n}\times {\mathbb{R}}^{r}\times \mathbb{R}\to \mathbb{R}$, $\boldsymbol{f}:{\mathbb{R}}^{n}\times {\mathbb{R}}^{r}\times \mathbb{R}\to {\mathbb{R}}^{n}$ and $\boldsymbol{g}:{\mathbb{R}}^{n}\times {\mathbb{R}}^{r}\times \mathbb{R}\to {\mathbb{R}}^{s}$ are ${C}^{1}$ -continuous with respect to all their arguments. We call $\left(\boldsymbol{x}\left(t\right), \boldsymbol{u}\left(t\right)\right)$ an admissible pair if $\boldsymbol{u}\left(t\right)$ is any piecewise control and $\boldsymbol{x}\left(t\right)$ is ${C}^{1}$ -continuous such that (A.5)–(A.10) are satisfied.

(Filippov-Cesari's Theorem, ^[60]) Suppose that there exists an admissible pair $\left(\boldsymbol{x}\left(t\right), \boldsymbol{u}\left(t\right)\right)$ and further that

(1) $\mathcal{U}$ is closed.

(2) $G\left(\boldsymbol{x}, t\right) = \left\{ {\boldsymbol{\tilde y}}\equiv \left(\boldsymbol{y}, {\boldsymbol{y}}_{n+1}\right):\boldsymbol{y} = \boldsymbol{f}\left(\boldsymbol{x}, \boldsymbol{u}, t\right), {\boldsymbol{y}}_{n+1}\ge F\left(\boldsymbol{x}, \boldsymbol{u}, t\right), \boldsymbol{g}\left(\boldsymbol{x}, \boldsymbol{u}, t\right)\ge 0, \boldsymbol{u}\left(t\right)\in \mathcal{U}\right\}$ is convex for all $\left(\boldsymbol{x}, t\right)\in {\mathbb{R}}^{n}\times \left[{t}_{0}, {t}_{1}\right]$.

(3) There exists a number $\epsilon > 0$ such that $‖\boldsymbol{x}‖ < \epsilon$ for all admissible pairs $\left(\boldsymbol{x}\left(t\right), \boldsymbol{u}\left(t\right)\right)$ and all $t\in \left[{t}_{0}, {t}_{1}\right]$.

(4) There exists an open ball $\mathcal{B}\left({\bf{0}}, \xi \right)\subset {\mathbb{R}}^{r}$ which contains the set $\mathcal{W}\left(\boldsymbol{x}, t\right) = \left\{\boldsymbol{u}\left(t\right)\in \mathcal{U}:\boldsymbol{g}\left(\boldsymbol{x}, \boldsymbol{u}, t\right)\ge {\bf{0}} \right\}$ for all $\boldsymbol{x}\in \mathcal{B}\left({\bf{0}}, \xi \right)$.

Then there exists an optimal pair $\left({\boldsymbol{x}}^{*}\left(t\right), {\boldsymbol{u}}^{*}\left(t\right)\right)$ to (A.4)–(A.10) with ${\boldsymbol{u}}^{*}\left(t\right)$ measurable.

For (A.1)–(A.3), the Filippov-Caseri's Theorem can be easily verified. Hence, the existence and uniqueness of Problem 2 is guaranteed.

Acknowledgments

The authors are grateful for the National Key Research and Development Plan (2019YFB1706502); the Natural Science Foundation of Shandong Province (ZR2020QG055); the National Natural Science Foundation of China (12102077); the Fundamental Research Funds for the Central Universities (DUT20YG125).

Conflict of interest

The authors declare no conflicts of interest.