Optimal vaccination strategies using a distributed model applied to COVID-19.

Optimal distribution of vaccines to achieve high population immunity levels is a desirable aim in infectious disease epidemiology. A distributed optimal control epidemiological model that accounts for vaccination was developed and applied to the case of the COVID-19 pandemic. The model proposed here is nonstandard and takes into account the heterogeneity of the infected sub-population with respect to the time since infection, which is essential in the case of COVID-19. Based on the epidemiological characteristics of COVID-19 we analyze several vaccination scenarios and an optimal vaccination policy. In particular we consider random vaccination over the whole population and the prioritization of age groups such as the elderly and compare the effects with the optimal solution. Numerical results of the model show that random vaccination is efficient in reducing the overall number of infected individuals. Prioritization of the elderly leads to lower mortality though. The optimal strategy in terms of total deaths is early prioritization of those groups having the highest contact rates.

Introduction

The control of the globally fast spreading severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) - the virus that causes coronavirus disease 2019 (COVID-19) - was more difficult than expected. Symptom based screening alone could not bring the pandemic under control. Accumulating evidence showed that SARS-CoV-2 is transmitted from persons without symptoms (Johansson et al. 2021) making infection control very difficult. In the absence of vaccines and effective drugs non-pharmaceutical interventions were the first line of defense against the pandemic. Population and individual based interventions such as isolation of infected individuals and quarantine of suspected infections were implemented widely. These measures led to a slowing down of the pandemic in the short-term. However, lifting of the measures led to resurgence of infections among people.

The remarkable speed of vaccine development became a powerful pharmaceutical tool for the control of the COVID-19 pandemic. Fast and optimal deployment of vaccination is expected to lead to a large fraction of the population with immunity, which will add to post-infection natural immunity of another part of the population. These two parallel processes will allow for reaching levels of immunity close to herd immunity, which will reduce future fatalities and severe cases. Vaccination aims at protecting individuals against infection and disease severity and reducing transmission of infected vaccinated individuals to others. The epidemiology of COVID-19 showed that severity of disease is variable among individuals, but increases with age (Davis et al. 2020). The age dependence was also observed for lethality with a strong concentration in the elderly population (Onder et al. 2020). On the other hand, the number of contacts also changes with age (Bhattacharaya et al. 2016), which influences the transmission of the pathogen. These characteristics point to the necessity of vaccine prioritization. This can be vaccination of those at highest risk of severe disease and potentially death. It can also be vaccination of those who transmit the virus most. Therefore, a fundamental question raised was how to optimize the vaccine allocation to maximize public health benefit. For instance, is it better to vaccinate younger adults or older people first? (Bubar et al. 2021) showed in a modeling approach that vaccination of younger adults may prevent the highest incidence whereas vaccination of elderly may reduce mortality. A mathematical modeling study explored different strategies for the prioritization of SARS-CoV-2 vaccines (Jentsch et al. 2021). Strategies that interrupt transmission such as vaccination of the younger population groups or vaccinating uniformly by age are more effective over time as herd immunity builds. They point to the role of timing of vaccination start. Vaccination of the elderly first was beneficial in the early stage of the pandemic. Late vaccination start would favor the employment of transmission interruption strategies over prioritizing the oldest in terms of preventing deaths. Vaccine uptake and effectiveness play also a critical role in a successful pandemic control including risks associated with early relaxation of non-pharmaceutical interventions (Moore et al. 2021). Looking at vaccine effectiveness scenario modeling showed that a vaccine with more than 50% effectiveness will lead to pandemic control, provided that a high percentage of the population is optimally vaccinated (Matrajt et al. 2021). According to this model the optimal solution to minimize deaths is to vaccinate the older population first, if vaccination effectiveness is low. For high vaccine effectiveness, transmission interruption by vaccinating high transmission age groups (usually younger people) was the optimal strategy.

In this work we adapt a previously developed epidemiological model (Kovacevic et al. 2022) to account for vaccination in order to implement several vaccination strategies and identify optimal solutions. The original distributed optimal control epidemiological model is based on a single integral equation having the number of new infections over time as the only dynamic variable. The basic approach aims at modeling the course of infection. Therefore, different subsequent sub-population compartments are considered, such as asymptomatic infected individuals, infected individuals with clinical symptoms (both isolated or not), recovered and dead individuals. The time since infection of individuals in any compartment is explicitly involved in the model, which therefore, becomes distributed and has the form of a system of integral equations. In Sect. 2 we present the fundamental characteristics of the original model and extend it to incorporate vaccination and describe vaccination scenarios. First, we consider random vaccination where susceptible individuals are randomly vaccinated with a certain rate in a certain period of time. We derive the total size of the vaccinated population and intensity before and at time t. We also extend the model with vaccination to account for heterogeneity of vaccination prioritized groups, we discretize the population of susceptible and recovered individuals and formulate a vaccination control policy. Minimizing the number of deaths is formulated as the control objective and the analysis of the optimization problem and is presented in Sect. 3. In Sect. 4 numerical results from plausible scenarios for COVID-19 control through vaccination are presented. These include, prioritization according to the severity of disease and risk of death that is prioritization of the elderly, transmission interruption by prioritizing younger adults, and comparison with the corresponding optimal solutions. The numerical results and the flexibility of the modified model in addressing multiple scenarios are discussed in Sect. 5.

Description of the model

In this section we introduce the model including vaccination policies as decision variables. In that we build on the setup in Kovacevic et al. (2022), which is modified to deal with several sub-populations consisting of people having similar levels of risk, say, due to age differences, and for which vaccination could be prioritized. Therefore, we start with a short review of the original model and then introduce the necessary modifications in detail.

The basic setup without vaccination

Kovacevic et al. (2022) split the population into compartments representing the stages of the disease that an individual may undergo during an epidemic:While transitions between the compartments are usually modeled by ordinary differential equations and transition rates, we use integral equations which allow to model explicitly the course of infection. Therefore, the mathematical formulation is based on the probabilities that an infected individual belongs to one or another compartment at a given infection age (time since infection).

S – susceptible

– asymptomatic, non-isolated

– asymptomatic, isolated

– symptomatic, non-isolated

– symptomatic, isolated

R – recovered

D – dead

These probabilities have been analyzed in the medical literature for many infectious diseases, whereas health authorities gain experience regarding the timing of quarantine and isolation measures.

The relative sizes of the compartments of infected individuals strongly depend on the infection age, as well as some key parameters such as the infectiousness and the contact rates. For this reason, the infection age plays an essential role in the model. We denote it by . Thus an individual infected at time will have infection age at time .

Kovacevic et al. (2022) assume that the environment is stationary (i.e. meteorological or behavioral changes are ignored1). Moreover, it is assumed there that recovered individuals remain immune after recovery. Both assumptions are kept throughout this work.

As in Kovacevic et al. (2022), we consider a finite time horizon [0, T], while the infection may have started earlier than time 0. It is assumed that all infected individuals either recover or die till infection age (say, 40 days). The following data related to the progress of the disease along the infection age are required. Each of them represents the fraction of all individuals infected at the same time, having a given status () at infection age :In accordance with the meaning of , we extend , , and , for .

– fraction of asymptomatic non-isolated individuals (

– fraction of asymptomatic isolated individuals ()

– fraction of symptomatic non-isolated individuals ()

– fraction of symptomatic isolated individuals ()

– fraction of recovered individuals (R)

– fraction of dead individuals (D)

Clearly, it must hold that for all .The above data allow to express the sizes of each of the compartments S, at any time by means of a single function of time, , which gives the number of new infections at any given time . Somewhat overloading the notation, denotes the size of the compartment at time () and similarly for S(t), R(t) and D(t). Moreover, for are considered as known initial data, therefore they appear separately in the right-most expressions below. We have (using the extensions of the fractions introduced above for ) thatwhere is the size of the susceptible population at time , is the size of the compartment of individuals who have been infected before time and have recovered till time , is similar but for the dead individuals.

Thus, the main unknown variable in the model is y(t) – the quantity of individuals who get infected (for the first time) at time t. According to the definition of , all these individuals either recover or die at most days after infection.

More data is needed for modeling the dynamics of y(t), some of which are time-dependent in order to incorporate possible distancing policy as in Kovacevic et al. (2022), namelyIn the derivation of the equation for y below, natural deaths and births are ignored, which is plausible if the duration of the epidemic is not too long or the natural demographic change is slow.

c(t)    – contact rate of susceptible and recovered individuals,

– contact rate of individuals from compartment , ,

– infectiousness of individuals from compartment , , of infection age .

It is assumed that the amount of newly infected individuals at time t is proportional to the number of contacting susceptibles, S(t), regarding their contact rate, c(t), and to the infectiousness per contact of the environment in which the contacts take place:The infectiousness of the contact environment is determined on the assumption of weighted random mixing. For that we introduce the cumulative infectiousness of each of the four compartments of infected individuals:Then takes the formwhereare the cumulative infectiousness of the contact environment and the total number of contacts per unit of time, respectively. Substituting the obtained expressions for , P and Q in (2.5) and using (2.1), (2.2) and (2.6), we obtain the equationwhereNotice that is the only unknown variable in the model and (2.8) is an evolutionary integral equation for its dynamic. All other functions that directly or indirectly appear in (2.8) are data that are assumed to be known.

We mention that the model can be further extended, say by including hospitalized individuals, hospitalized in intensive care, etc., provided that data similar to are available for these compartments. Such extensions, however, do not change the structural form of the transmission dynamics (2.8).

Modeling random vaccination

Using the basic notations from the previous subsection we extend the model step by step. First, we consider the simplest vaccination scenario where susceptible individuals are randomly vaccinated with a certain intensity and in a certain period of time. The vaccinated individuals are treated in the same way as recovered ones. It is assumed that the immunity after vaccination, as well as after recovery, is permanent (considering waning immunity will be subject of future work). Also, if there is a delay between vaccination and immunity or several doses have to be given over time, we model the effective vaccination as occurring at the time when immunity is reached.

Vaccination begins at time . Since we assume permanent immunity, the vaccination becomes redundant after a while. We assume that it ends at time , and that in the period of time there are enough vaccines, vaccination facilities, and individuals willing to be vaccinated, so that one can achieve a given vaccination intensity v(t) at any time . This means that the total size of the vaccinated population prior to time isThen regarding the vaccination, the basic model remains essentially the same, with the following modifications: the expressions for S(t) and R(t) in (2.2) and (2.3) are replaced withand this expression for S(t) has to be used in (2.5). The multiplier does not change, because P(t) depends on only, and the terms w(t) cancel in the expression for Q(t). Since the vaccinated individuals are treated as recovered, the basic Eq. (2.8) remains the same, with the only difference that has to be replaced with , whereIn this modification of the model it is implicitly assumed that the vaccinated individuals are chosen “randomly”, that is, disregarding the possible heterogeneity of the population with respect to disease-related mortality and contact rates. An enhancement of the model, taking into account the heterogeneity of the population follows in the next subsection.

Extension of the vaccination model involving population heterogeneity

In order to enable prioritized vaccination of certain population groups we split the population compartments S, R and in the basic model into m groups:The affiliation of a susceptible individual to a group is preserved after infection, recovery or vaccination. Each of the groups has specific parameters , , , , , , . For example, is the fraction of asymptomatic non-isolated individuals of infection age in the group j, among all the individuals in the same group and infection age , is the contact rate of non-isolated infected individuals of type j during the asymptomatic period, etc. Clearly, the identitieshave to be fulfilled. We also split the new cases into m groups , , where , are known initial data for each group. The meaning of and is clear from the notations. The group-specific form of expression (2.1) is straightforward:Following the group splitting, we denote by the intensity of vaccination at time t of members of group j. Correspondingly, denotes the size of the vaccinated till time t subpopulation of group j. Then the formulae (2.6) and (2.10) take the group-wise formAssuming random mixing of all groups, we have similar to (2.5) and (2.7) thatwhere P(t) and Q(t) have the same meaning as in the basic model, but are now given bySubstituting the above expressions in (2.14) we obtain the systemwhereA vaccination policy is any allocation of the available vaccination capacity, , to the groups of susceptibles. That is, any m-tuple of measurable functions satisfying the constraintsis considered as feasible control (vaccination) policy. Thenis the number of vaccinated before time t individuals from the j-th group, which appears in (2.15).

As mentioned above, Eq. (2.15) is based on the assumption of random mixing over all (pooled) groups. This approach can be easily enhanced by considering mixing with various contact rates between and within groups. To this end, we consider contact rates for the susceptible (recovered) individuals in group i with individuals of group j. In addition, contact rates for the infected compartments are needed. In the simplest case, they can be calculated as , where is a group specific factor, modeling the reduction of contacts.

Such an approach can be used to model the relations between different behavioral groups in more detail, accounting for the fact that the depth of relations between varying pairs of groups may substantially differ. Elderly people in nursing homes or children are examples of groups with intense contacts within the group and much lower outside the group.

The equations for the newly infected within groups becomewhereThe main drawback of the extension involving various contact rates between and within groups, is that it requires more data that are hardly available. Therefore we will not use it in the present paper.

Optimal vaccination policy

As explained in the preceding section, the vaccination intensities, , will be used as control policies. Feasible controls (allocation policies for the available vaccines) are the measurable functions satisfying the constraints (2.17). The controls influence the dynamics (2.15) through the cumulative quantities in (2.18).

Various control objectives may be reasonable. One goal of the vaccination control may be to ensure (together with additional policies such as social distancing) that the number of individuals needing hospitalization does not exceed the capacity of the hospitals. This constraint takes the formwhere is the fraction of symptomatic isolated individuals of infection age who need hospitalization, and H(t) is the capacity of hospitals at time t. Clearly, the first term in the above expression is just a known number.

The constraint (3.1) may be regarded in the context of optimization involving contact restrictions, where the economic consequences of the latter should be minimized. The economic consequences of “lock-down” have been recently investigated in several papers, e.g. (Acemoglu et al. 2000; Bloom et al. 2000; Kovacevic et al. 2022; Caulkins et al. 2021). However, in the present paper we focus on the purely health benefits of vaccination, therefore we set the minimization of deaths as objective. Formally, the problem reads asWe remind that is the fraction of dead individuals of infection age in the j-th group, among all infected individuals of the same infection age and group. Since the first term in the above expression is independent of the control, we can reformulate the problem assubject to the state dynamics (2.15) and the constraints (2.17)–(2.18).

For convenience we reformulate the problem in a more general way as follows. The state is a function , the control is a measurable function . We also define an aggregated state and control . The controlled dynamics is described by the equationswhere , . In addition, there are control constraints of the formwhere .2 The objective functional to be minimized will have the formwhere is continuous and denotes the scalar product in . Thus we encounter an optimal control problem for a specific class of integro-differential systems.

Further we denote by the set of all measurable functions (feasible controls) v satisfying (3.6), and by the set of all functions w resulting from (3.5) for some . In addition, denote by W the set of all values that w(t) may take when .

Obviously, our epidemiological problem is a special case of problem (3.3)–(3.7). Indeed, set . For convenience we denote the components of z by . Then we recast the function in (2.15) aswhere z is given by (3.4) with the following matrix-function P (skipping the arguments in (2.16))The objective integrand iswhere means transposition.

In order to facilitate the analysis of problem (3.3)–(3.7), in the next paragraphs we introduce a few notations and make some assumptions.

The notations and will be used for the space of all measurable (vector-) functions on [0, T] that are square integrable, respectively essentially bounded. The norm in will be denoted by , while a weighted norm in will be used in the proof of the next proposition 3.1.

For any , let be a map on with values in the space of all measurable functions on [0, T] defined asWe assume that there exists a closed bounded convex subset of which is invariant with respect to any , meaning that for every and . Moreover, we assume that the function f(t, z, w) is continuous in t and differentiable in (z, w) with Lipschitz continuous derivatives (uniformly with respect to ) in a neighborhood of the set , whereBefore we begin with the analysis of problem (3.3)–(3.7), we comment on the verification of the assumptions made above for the particular optimal vaccination problem (3.2), (2.15), (2.17)–(2.18). There are two problems in the verification of the assumptions concerning the function f in (3.8). One is to ensure that the second multiplier in the numerator remains positive (that is, no group totally dies out), the second is that the denominator does not approach zero (that is, the amount of contacting individuals stays above certain positive level). Both requirements are completely plausible from epidemiological point of view, but require formal assumptions that we skip in this paper; essentially they are similar to those formulated in Sect. 2.2 of Kovacevic et al. (2022).

Now we return to the general problem (3.3)–(3.7).

Proposition 3.1

On the assumptions made, for every feasible control the system (3.3)–(3.5) has a unique solution .

Proof

As usual, the proof uses the contraction mapping theorem; we present it in detail because system (3.3)–(3.4) does not fit to the standard format Volterra integral equations.

In the space we consider the weighted normFor an arbitrarily fixed and the corresponding w in (3.5) we shall prove that the map is contractive with respect to this norm, provided that the number is chosen sufficiently large. This implies the existence and uniqueness claims of the proposition.

Let L be the Lipschitz constant of f with respect to z in the set Z. For any we haveClearly, by choosing sufficiently large we may ensure that the operator is contractive on with respect to the norm in . Thus it has a unique fixed point , and it (together with z given by (3.4)) is the unique solution of system (3.3)–(3.4) in for the fixed v.

Proposition 3.2

The optimal control problem (3.3)–(3.7) has a solution.

If is a minimizing sequence with and then, due to the weak compactness of and , one may extract a subsequence which is weakly convergent in to some . Then due to (3.4) and (3.5), the sequences and converge point-wise to and , respectively. Then converges point-wise to f(t, z(t), w(t)). This implies that solves (3.3)–(3.5) and must be an optimal solution, because one can pass to the limit in (3.7).

It turns out that the functional F in (3.7) is Fréchet differentiable in , which is to be expected due to the linear structure of the problem with respect to v and y. Details of the derivation of the next result are presented in Kovacevic et al. (2022), Appendix.

Proposition 3.3

On the assumptions made above, the objective functional F(v) is Fréchet differentiable in at every , and its derivative has the representationwhere w and z correspond to v, and is the unique solution of the Volterra integral equation

Now, let with and be an optimal solution of problem (3.3)–(3.7). Denote by V the set of feasible control values, that is . As a consequence of Proposition 3.3 we obtain the following necessary optimality condition.

Corollary 3.4

If is an optimal solution of problem (3.3)–(3.7), then there exists a unique solution of (3.9), with (z, w) replaced with , such thatwhere is the usual normal cone to the convex set V at the point v.

We mention that the existence of the derivative of the objective functional with respect to the control and the constructive way of its calculation open the door for efficient numerical algorithms, such as the gradient projection method in the control space. However, in the present paper we do no discuss the numerical aspects of the problem.

Let us return to the epidemiological problem (3.2) subject to the state dynamics (2.15) and the constraints (2.17)–(2.18). Since is the first order increment of newly infected individuals resulting from an increment of the number of vaccinated individuals, the expression for in Proposition 3.3 means that is the marginal “cost” (increment of dead people) resulting from a unit increase of the new infections at time t. In this sense, is the “shadow price” of new infections at time t at the optimal solution.

Notice that the set of feasible control values, V, is a simplex, at any vertex of which at most one of the components may be non-zero and equals . Since the variational inequality (3.10) is linear (meaning that the left-hand side is independent of v), one can expect that the optimal control has a switching structure: at any time only members of one group are vaccinated with the maximal intensity . Of course, presence of singular arcs (where the optimality condition does not uniquely define a control) cannot be excluded a priori. However, in our experiments with the COVID-19 model we always encounter purely switching structure of the optimal control (see next section).

Remark 3.5

As mentioned above, the objective functional (3.2) has a purely “humanitarian” meaning. The cost of the vaccination is ignored, as well as the economic benefits from the reduction of the diseased population due to vaccination. Including a strongly convex cost of vaccination with a “small” weighting coefficient can be viewed as Tikhonov regularization of the problem. The optimal vaccination policy in this case will no longer be of bang-bang type, however, it will approximate in the space the optimal policy for the case of zero vaccination cost.

Case studies

We use the presented COVID-19 vaccination model to compare the optimal vaccination strategy for age groups with several plausible vaccination strategies. The model is parametrized based on available epidemiological and sociological data as described below, see also (Kovacevic et al. 2022) for a deeper discussion of the course of infection.

Basic parametrization

In our vaccination model we consider five age groups of different individuals: years, years, years, years and over 80 years. The considered age structure depicts one young age group, two groups of working people and two groups of retired people. This allows us to use groups specific sizes, contact and mortality rates.

The size of these groups is derived from official information on five years groups at Statistics Austria (2020). The different contact rates of our groups were modeled based on Bhattacharaya et al. (2016) Fig. 1, depicting the average number of alters for egos of varying age. For each of the considered age groups, the proportion of cumulated alters (compared with cumulated alters over 100 years) is used as a modifying factor for the overall contact rate c. Here it is assumed that the contact rate is roughly proportional to the age specific number of alters.

Fig. 1

Comparison between the percentage of infected individuals (left) and comparison between mortality cases (right) with different vaccination strategies

For characterizing the mortality, we use (Signorelli and Odone 2020). Note that the age groups in the mentioned paper are different from the groups in our case study. Therefore, the mortality is recalculated by accounting for the size of five years subgroups within our subgroup population. Here, it has to be assumed that mortality is constant over all five years groups within any age group.

In our model only individuals with symptoms may die. In order to describe the morality depending on time since infection , we use the mortality function estimated by Verity et al. (2020). This raw mortality function is normalized such that the area below the function is equal to one. The functions then can be constructed by multiplying the normalized mortality function by the proportion of eventually symptomatic individuals within the class of infected in group j and with the proportion of fatally ill individuals within the class of eventually symptomatic individuals in group j.

In Table 1 we provide the values describing the group specific parameters. Note that the contact rates in Table 1 are used as factors: they are used for deriving the group contact rate by multiplication with the overall contact rate c of the population.

Table 1

Parameters for age specific groups

Groups	0–18	18–30	30–65	65–80	80+
Fraction of population	0.1933	0.1257	0.4907	0.1370	0.0532
Contact rate (factor)	0.5271	1.6209	1.2250	0.7982	0.7011
Case fatality rate	0.0010	0.0010	0.0369	0.2232	0.3460

The class of infected individuals at given time t consists of individuals who show symptoms over some time and of completely asymptomatic infected individuals, who recover without having shown any symptoms.

Finally, the available epidemiological information in the course of the infection can be used to construct parameter functions and , similar to the parameter functions and c,  described in Kovacevic et al. (2022).

In Table 2 we provide values of the parameters that are used for the numerical experiments and are different of the previously used in Kovacevic et al. (2022).

Table 2

Parameter values

Overall probability of infection per contact	0.1429
Overall contact rate	12.978
Reduction of contact rate for isolated individuals	0.09
Reduction of contact rate for non isolated symptomatic individuals	0.97
Infected without symptoms	0.425
Percentage of isolated asymptomatic cases	0.60
Percentage of isolated symptomatic cases	0.85

Numerical results

Below we present numerical results showing the evolution of the epidemic without and with implementation of different vaccination policies. The vaccines at time t are only administered to individuals that do not show symptoms, i.e. individuals of the different age groups who do not belong to classes and .

We consider the following vaccination scenarios.

Case 1: no vaccination is applied during the duration of the disease.

Case 2: random vaccination is administered to the population with an exception of the age group .

Case 3: the vaccination begins with the vaccination of the elder population, that is, the two age groups of and . After all suitable for vaccination individuals of these groups are vaccinated, the vaccination continues by applying the obtained by numerical optimization vaccine policy to the rest of the groups.

Case 4: the numerical calculated optimal vaccination strategy is applied.

Additional experiments regarding the optimal vaccination order in our model (Case 4) are also presented, exploring the importance of the contact rates of the groups.

For all scenarios we consider a time horizon from 0 to 455 days. The time interval on which the vaccination control v(t) is applied is [150, 250].One can easily extend or change the vaccination administration interval .

The vaccination period is chosen to contain the majority of the infected cases in our experiments and also to demonstrate the differences in the vaccination strategies. Having this period shifted more to the left of the horizon, i.e. vaccination is administered earlier, has an even greater effect of reducing infections and mortality. This is due to the fact that vaccinated people are protected against future infections.

In the experiments the second constraint in (2.17) reads as , i.e. the daily vaccination capability is of the population. The value is chosen to represent somewhat realistic constraint on the daily capability of administering vaccines either from supply, personnel considerations or people’s willingness to be vaccinated.

We compare the results from experiments with different case studies. On Tables 3, 4 and 5 are shown the vaccination policies for the different vaccination approaches. During each vaccination period the maximum daily vaccination capacity is administered to the corresponding age group. In the case of random vaccination, Table 4, the vaccination is administered simultaneously to all the groups except the youngest age group, that is why the groups vaccination period coincides with the period.

Table 3

Control policy for vaccination of first

Age group	Vaccination period	% of vaccinated within the group	% of vaccinated
65+	150–213	100%	16.9%
30–65	214–250	39.68%	12.9%

Table 4

Control policy for random vaccination

Age group	Vaccination period	% of vaccinated within the group	% of vaccinated
18–30	150–250	47.7%	4.2%
30–65	150–250	47.7%	18%
65–80	150–250	47.7%	5.5%
80+	150–250	47.7%	2.1%

Table 5

Control policy for optimal vaccination

Age group	Vaccination period	% of vaccinated within the group	% of vaccinated
18–30	150–190	100%	10.9%
80+	191–210	100%	4.9%
65–80	211–240	100%	11.2%
30–65	241–250	7.3%	2.8%

In Table 3 we see that after 63 days of vaccination all elderly people are vaccinated. After that our optimal solution suggests to start vaccinating the age group of years.

Different is the case of the optimal vaccination policy shown in Table 5. The obtained optimal solution suggests to begin the vaccination with the highest contact rate group of individuals years. After the group of young working is fully vaccinated on day 190 the solution switches policy to the people with highest mortality rate () and again after this group is vaccinated the vaccination continues with the group of years. The vaccination ends with the vaccination of the rest of the working population, years. We explain that the prioritization of the young working group at the beginning of the vaccination interval is due to the high impact of the contact rate parameter in our model, i.e. although the elderly population has the highest mortality rates, we try to minimize the total number of fatality cases among all groups.

In addition, in Tables 6 and 7, we show optimal solutions obtained when changes are made to the contact rates of group years and group years. For the experiment in Table 6 the contact rate factor for group is reduced to the value of group years, i.e. from 1.6209 to 1.2250. Here we can observe that the order of vaccination is changed, the vaccination begins with the group years. In contrast to the previous observed optimal behavior, the vaccination ends before all individuals in the group are vaccinated. This is due to the size of this group which is significantly larger than the rest of the groups. Afterwards the strategy again continues with the vaccination of the two elderly groups.

Table 6

Control policy for optimal vaccination with reduced 18-30 contact rate

Age group	Vaccination period	% of vaccinated within the group	% of vaccinated
30–65	150–210	52.3%	16.6%
80+	211–230	100%	5.2%
65–80	231–250	71.2%	8.1%

Table 7

Control policy for optimal vaccination with increased 80+ contact rate

Age group	Vaccination period	% of vaccinated within the group	% of vaccinated
80+	150–168	100%	4.7%
65–80	169–208	100%	11.3%
30–65	209–250	43.7%	13.9%

For the experiment in Table 7 the contact rate factor for group is increased again to the value of group years, i.e. from 0.7011 to 1.2250. For this example we observe again changes in the order of vaccination. This time the vaccination begins with group , and — after all group members are vaccinated — continues with group . Finally it ends with the vaccination of group .

It should be noted that varying one parameter changes the overall disease evolution over time. The obtained optimal vaccination solutions with respect to the various contact rate factors are compared only with the changed optimal vaccination order.

Remark 4.1

Shrinking the vaccination period plays also a role for the optimal vaccination. Having a small interval for vaccination prioritizes the most vulnerable, i.e. the elderly population, to be vaccinated first. In our case increasing the end of the vaccination horizon preserves the obtained strategy to vaccinate first the group with the highest contact rate.

On Fig. 1 are represented the results for the different vaccination strategies showing the total percentage of infected individuals and the total percentage of mortality cases for all groups. From the left plot of Fig. 1 one can observe that the random vaccination strategy is more efficient in reducing the overall infected percentage compared to strategy of vaccination the elderly population of people. Although, the case fatality rates of age groups and is significantly higher than the other groups, we can observe that by prioritizing the vaccination of the the total mortality is comparable to the random vaccination. In the case of the optimal vaccination policy (Table 5) on Fig. 1 it’s shown that this strategy is superior to the other cases, both in reducing the number of infected individuals and reducing the overall fatality cases.

During the vaccination period the same number of vaccines are administered to the population with different results. Figure 2 shows the evolution of the population compartments during the evolution of the disease. One can observe the changes of the levels not only for the infected individuals but also for the recovered. After the vaccination period [150, 250], between and of the population is vaccinated or recovered from the disease depending on the strategy. The rest of the population consist mainly of susceptible individuals, but small portion of infected individuals (that may show no symptoms) remains. This means one can continue to vaccinate even after the considered vaccination period.

Fig. 2

Comparison between the different vaccination strategies with respect to the population compartments

Discussion

In this study, we build on a previously developed epidemiological model by including vaccination strategies. Several vaccination options are investigated, such as random vaccination or prioritization of population groups according to their age classification. The latter is of particular interest for the application of the model to COVID-19. It is well known that age is a critical population vulnerability factor. The risk of severe illness with COVID-19 increases with age where the elderly population is at highest risk (Wu et al. 2020). In an optimal control framework we develop vaccination policies, that is, the allocation of the available vaccination capacities to the susceptibles to be vaccinated. Control objectives were chosen here based on the most desirable effects such as to keep the number of hospitalizations such as it does not exceed the capacity of the hospital and therefore does not overwhelm the health care system. Another control goal may be to minimize the number of fatalities. This was also the control objective in the COVID-19 case study.

We formulate the model for random vaccination with heterogeneous populations. Other options such as mixing with different contact rates between and within population groups such as elderly, supporting nurses, physicians and the general population could be incorporated in the model in a straight forward way. While any relaxation of the assumptions and expansion of the model is mathematically doable, the availability of corresponding data is rather difficult. The epidemiological optimal control policy is derived from an optimal control problem for a specific class of of integro-differential equations system with certain characteristics such as control-affinity. Its analysis indicates practical decisions such as vaccination of one group at a time.

The numerical results are based on plausible vaccination strategies for the case of COVID-19. The random vaccination is applied to the population with the exception of age group 0-18 since the latter group is by far the least vulnerable to develop severe disease although it contributes to transmission. In another scenario we prioritize vaccination of the elderly population, a strategy that has been widely followed. Finally, we calculated the optimal vaccination strategy. In comparison to an non-vaccination scenario where the epidemic runs without intervention, vaccination in whatever scenario had a clear positive effect in reducing the number of cases and shortening the duration of the outbreak. In the scenario of prioritization of the elderly the model results show that the optimal solution continues with the age group . Interestingly, in the optimal vaccination policy vaccination begins with the population groups having the highest contact rates. This may have implications when designing the vaccination strategy. If other restrictions such as ethical reasons do not speak against changing the order to an optimal scenario. Thus, a clinically advisable scenario in not necessarily also epidemiologically optimal. Clearly, changes of the contact rates may lead to changes in the optimal order at which different age groups should be vaccinated.

Although the presented model can be easily extended in several directions (e.g. by including additional compartments, medication, imperfect effect of vaccination, etc.), there are substantial limitations to its applicability. Of course, as for all complex epidemiological models, a main problem is the data availability. In addition, we assume permanent immunity, which may be plausible only on the short run. Consideration of waning immunity (after vaccination or disease) is a subject of a paper in preparation by the same authors. In addition, we mention that the control strategies considered in this paper, including the optimal ones, may not be fully practical. Such difficulties arise in particular with frequently changing policies or policies requiring vaccination of all members of a group. Nevertheless, the information provided by such strategies gives a hint for what policy should be targeted by the health care authorities.