Dynamical intervention planning against COVID-19-like epidemics

COVID-19 has got us to face a new situation where, for the lack of ready-to-use vaccines, it is necessary to support vaccination with complex non-pharmaceutical strategies. In this paper, we provide a novel Mixed Integer Nonlinear Programming formulation for fine-grained optimal intervention planning (i.e., at the level of the single day) against newborn epidemics like COVID-19, where a modified SIR model accounting for heterogeneous population classes, social distancing and several types of vaccines (each with its efficacy and delayed effects), allows us to plan an optimal mixed strategy (both pharmaceutical and non-pharmaceutical) that takes into account both the vaccine availability in limited batches at selected time instants and the need for second doses while keeping hospitalizations and intensive care occupancy below a threshold and requiring that new infections die out at the end of the planning horizon. In order to show the effectiveness of the proposed formulation, we analyze a case study for Italy with realistic parameters.


Introduction
The ongoing COVID-19 pandemics, with its huge toll in terms of deaths and economic damage, represents an unparalleled global threat to human society as a whole. As of May 2021, reportedly more than 155 million COVID-19 cases have been identified, with more than three million deaths [1]. To face such a threat, governments have initially reacted via non-pharmaceutical interventions, i.e., by enforcing strict social distancing [2][3][4][5]. Then, with an unprecedented effort by human society as a whole, a wide variety of vaccines have been developed in a remarkably narrow time span [6][7][8][9][10]; such a rapid development has been possible also due to the extensive reliance on bioinformatics [7] and artificial intelligence [11]. Notably, the effectiveness and geographical distribution of such a plethora of vaccines has proven to be highly heterogeneous (e.g., see [12] and references therein). Notice that, to date, no universally acknowledged cure has been identified; the case of the debate regarding therapies based on Remdesivir represents an illustrative example in this sense [13,14]. Therefore, to date, the control knobs available to governments amount to just social distancing (e.g., lockdowns, limiting affluence to shops, wearing masks, etc.) and vaccination. In the literature, optimization tools to face epidemics, such as OpenMalaria [15,16] and STDSIM [17], have proved their effectiveness in analyzing and planning illness containment [18][19][20]. However, in the case of COVID-19, the unavailability of reliable information with adequate level of detail requires reliance on simpler approaches. Among several other options, SIR models [21,22] represent a reasonable choice in terms of predictive capability and simplicity. Interestingly, such models have been applied to describe different epidemics such as SARS [23] Influenza A (H1N1) [24], Measles [25] and Hepatitis C [26]. Moreover, compartmental models in general have proven useful to model quite different epidemics scenarios e.g., Chlamydia trachomatis, antiviral treatment in the case of HIV, nosocomial infections and transmission of antibiotic-resistant pathogens [27]. For this reason, such model allows to understand the available degrees of freedom, i.e., the policies that can be put in place to react to the epidemics, even in the absence of detailed quantitative predictions [28]. Indeed, several control [29][30][31] and optimization [29,[32][33][34][35] approaches have been developed, based on simple epidemics models such as the SIR. In particular, it is worth mentioning: [35], where a coarse-grained and static optimization framework for selecting the amount of vaccines to be allocated to different population classes with the aim of ending the epidemics is given; [36], where the authors focus on vaccination of essential workers; [37] where the authors allocate vaccines to age classes in order to optimize cost functionals such as deaths or hospitalization, but do not guarantee the end of the epidemics at the end of the considered time horizon nor consider different vaccines with different efficacies at the same time. Other examples include [38,39], where optimization of the supply chain underlying the vaccine delivery is considered. Notably, based on epidemiological data from the UK together with estimates of vaccine efficacy, [40] provides a framework to conduct "what if" analysis of the evolution of the disease under different interventions. However, while being beneficial for high-level policy making, does not translate into an operative plan for the administration of pharmaceutical and non-pharmaceutical interventions. Moreover, the workin [41] provides a model predictive approach to the problem where interventions at time t are selected based on the foreseen effect in the next future, following a "receding horizon" approach. However, although quite detailed, the epidemics model underlying such framework amounts to a single population of individuals and does not distinguish between age classes.
To the best of our knowledge, to date no formulation is able to support the fine-grain, operative and dynamical planning of interventions, while accounting for a wide range of phenomena that are peculiar of epidemics like COVID-19, e.g., different types of vaccines, the need for a second dose, the capacity of the healthcare system in terms of regular and intensive care hospitalization, the availability of vaccines in batches at selected time instants.
In this paper, we fill this gap by providing a novel optimization formulation that aims at implementing an optimal intervention plan to fight newborn epidemics like COVID-19, i.e., epidemics characterized by two key factors: (i) high infection rate and (ii) high stress posed on the healthcare system and/or society in terms of intensive care occupancy, deaths or economic consequences. Specifically, we first develop a modified discrete-time SIR model for heterogeneous population classes (e.g., age or geographical classes) that accounts for the effect of social distancing and vaccination. In more detail, we assume the ratio at which individuals in two classes infect each other can be reduced by enforcing tailored social distancing measures. Moreover, we consider several types of vaccines, each characterized by their efficacy as well as the delay required for the vaccination to be effective. In this view, we assume that, after the vaccine takes effect, a fraction of the population becomes immune. Notice that we explicitly take into account the possibility to plan for a first and a second dose of the vaccines.
Based on the proposed variation of the SIR model, we develop an optimization formulation that aims at planning social distancing measures and vaccinations at the level of the single day in order to reach the end of the epidemics and to guarantee that the state reached is robust, in that new infections die out. In doing so, we enforce constraints accounting for aspects such as the need of a second dose, the delayed and partial effect of multiple types of vaccines, the requirement that congestions of the healthcare system are avoided.
Briefly, we show that the proposed formulation amounts to a Mixed Integer Nonlinear Programming (MINLP) problem. We point out that the term MINLP is used to denote a class of optimization problems where some of the variables being selected are constrained to be integer-valued while some other can be real-valued (i.e., Mixed Integer (MI)), the objective function and/or some of the constraints are nonlinear functions (i.e., Nonlinear (NL)). In the context of optimization, "program" or "programming" (P) can be regarded as a synonym of optimization problem or optimization formulation. Notably, the class of MINLP problems is known to be computationally hard to solve (e.g., see [42]); therefore, exact solution of the problem is a nontrivial task, thus calling for approximation strategies to be put in place.

Notation
We denote vectors by boldface lowercase letters and matrices with uppercase letters and we refer to the (i, j)-th entry of a matrix A by A ij . We represent by 0 n and 1 n vectors with n components, all equal to zero and to one, respectively. Moreover, we use 0 n×m , 1 n×m to denote an n×m matrix with just zero and one entries, respectively. We use square brackets to denote the arguments of a function, e.g., we use f[x, y] to denote a function with arguments x and y. For the sake of brevity, where understood, we abbreviate a function of one or more arguments by f [�]. Given a vector x 2 R n we denote by diag[x] the n × n diagonal matrix having diag[x] ii = x i , for all i 2 {1, . . ., n}. On the same footings, given a matrix A 2 R n � R n we denote by diag[A] the n dimensional vector having diag[A] i = A ii , for all i 2 {1, . . ., n}. We denote by � the Hadamard (i.e., entry-wise) matrix product between matrices A and B with the same dimensions, i.e., the matrix C = A � B is such that C ij = A ij B ij ; analogously, the Hadamard product between c = a � b between two vectors a,b is such that c i = a i b i . Remember that Hadamard product is commutative; moreover, notice that

SIR epidemics model
In this section we briefly review the SIR Epidemics model; the interested reader is referred to [35,43] for further details. Let us consider a population of N individuals divided in n classes (e.g., by age or geographical area); we denote by N ℓ the population in the ℓ-th class with where B is the n × n transmission matrix, B ij being the rate at which a susceptible individual of class i meets an infectious individual of class j and becomes infected, while the vector g 2 R n collects the rates γ i at which infectious individuals in the i-th class are removed from the Modeling assumptions and limits of the SIR model. The simplicity of the SIR model allows to design a scenario based on a limited number of parameters; it is thus one of the first models used to understand newborn epidemics. However, SIR models can sometime oversimplify the complex disease process. As an example, SIR models imply "full mixing", i.e., the assumption that all individuals in the population are equally likely to be in contact with each other. To this respect, the heterogeneous SIR corrects such an issue by introducing classes and considering the heterogeneity in their contact rate. Also, we have employed a simplified SIR model with fixed populations, although in the original formulation it could account also for migration, births or deaths; such an approach is justified in the initial phase of an epidemic, where the time horizon is limited and variation in population can be disregarded. When an epidemic becomes endemic, SIR models can be easily extended to SIRS models where recovered individuals can become again susceptible. Furthermore, if there is a waiting time for an infected person to become infectious, SIR models can be extended to SEIR models by introducing an extra compartment E (i.e. "exposed") that accounts for such an issue. In a "receding horizon" approach, where the model parameters are periodically adjusted to reflect the evolving knowledge on the epidemic, it is reasonable to resort to the SIR model during the first iteration, since its parameters are the simplest to estimate. Eventually, at later iterations, it is possible not only adjust the parameters, but also to switch to more sophisticated models (SEIR or even SEIRS if the situation becomes endemic) with the proceeding of time. Our framework easily allows to switch from SIR to SEIR or SEIRS model just by adding extra compartments; notice that, by defining by the fraction of infected (i.e., exposed or infectious) individuals, the constrains ensuring the dampening of the epidemic for all classes (i.e., @a[t]<0) retain the same simple linear form of Eq (2), i.e., they depend only on the fraction of infectious individuals [44].

Modeling interventions within the SIR model
In this section we modify the SIR model in order to explicitly account for possible interventions, namely, social distancing measures (e.g., adoption of personal protection equipment (PPE) and lockdowns) and vaccination. In view of later developments in the paper, it is convenient to first express the SIR model in discrete-time form. Notice that exact discretization of a nonlinear differential equation However, given the complexity of the above exact method, it is convenient to consider an approximated relation. In particular, in this paper we resort to the Euler forward approximation, i.e., we set In other words, we approximate the continuous-time SIR model in Eq (1) via the following discrete-time equations and we point out that, to avoid numerical instability as a result of the discretization, we choose Δt = 0.01[day] for the parameters used in our case study (Other approaches allowing larger step size without causing instability could be considered, such as Euler backward integration, trapezoidal integration or Runge-Kutta methods; however, in this paper we opted for the Euler forward integration for the sake of simplicity). Let us now incorporate two different types of intervention in the above discrete-time SIR model, accounting for the adoption of social distancing measures and for vaccination. Notably, in the following, we consider interventions such as vaccinations and social distancing measures that are planned at the level of the single day; as discussed next, such interventions will reflect in the discrete-time SIR model by assuming that the interventions remain constant over the day. As a consequence, such daily interventions will be indexed on a daily basis, while the variables in the discrete-time SIR model are indexed by hΔt. Moreover, we use the iterator k to denote the k-th day and, where needed, with a slight abuse of notation we use the iterator k to denote the value assumed by a variable of interest at the end of the k-th day, e.g., s[k] = s[hΔt], with h = k/Δt, while we point out that the day corresponding to the time instant hΔt is given by bhΔtc.
Social distancing interventions. In order to model the effect of social distancing measures in the SIR model, we observe that interventions such as lockdowns, limiting access to shops or imposing the adoption of PPEs, has the effect to reduce the rate at which susceptible individuals meet infectious individuals and/or become infected, modeled by the coefficients B ij . In order to model the effect of social distancing measures at the k-th day, let us define the social distancing intensity E ℓj [k] 2 [0, e], where e is the maximum allowed intensity and e � 1, as the intensity of the social distancing measures put in place for the ℓ-th and j-th class, e.g., or, in a compact form Vaccination. Let us now model the effect of vaccination on the discrete-time SIR model. In particular, we assume m different types of vaccines are available and we assume each vaccine j has an efficacy η j 2 [0, 1] after a single dose, while after the second dose the efficacy rises to Notice that the second dose is not required for all types of vaccines; we model this aspect by resorting to a coefficient if the second dose is required for the jÀ th vaccine 0; otherwise: ( Moreover, for each type j of vaccine we assume a time window of t I j days is required for the vaccine to take effect after the first dose, while we use w j � t I j to denote the time window between the first and the second dose (if required) and t II j to denote the time window between the administration of the second dose and the reach of complete effect. In other words, an administration of vaccine j on the k-th day has an initial effect on day k þ t I j and a complete effect on day k þ w j þ t II j , while τ I , τ II are the times estimated from pharmacological trials during which vaccinated individuals are still exposed to the infection. In order to model the effect of vaccination, let us indicate with X ℓj [k] the units of first doses of vaccines of the j-th type that are injected to the ℓ-th class of population at the k-th day and let X½k� 2 N n�m be the matrix with integer entries collecting such variables. Moreover, let Y ℓj [k] denote the amount of units of second dose of vaccines of the j-th type that are injected to the ℓ-th class of population at the k-th day.
In this view, the contribution Δr ℓ [k] at day k to the number of removed individuals belonging to the ℓ-th class as a result of vaccination satisfies i.e., the contribution of the j-th vaccine to Δr ℓ [k] corresponds to the fraction of individuals that were vaccinated t I j days before with the first dose of the j-th vaccine, weighted by its efficacy η j , plus the fraction of individuals that were vaccinated t II j days before with the second dose of the j-th vaccine, weighted by the residual efficacy Δη j . In other words, Δr ℓ [k] is given by where, for the sake of consistency, we assume X[�], Y[�] are zero when their argument is negative.
Finally, we assume that i.e., if required (ϕ j = 1), the units of second dose of the j-th type of vaccine injected at the k-th day must not trespass the units of first dose injected χ j days before; otherwise (ϕ j = 0), no second dose is considered. Notice that, being Eq (8) an inequality, the second dose is not mandatory, and it is possible to implement policies where only a fraction of the population receiving the first dose receives also the second as suggested by the UK study SIREN [45].
Resulting SIR model. To conclude the section, let us show the expression of the discretetime SIR model where the above interventions are explicitly considered. In particular, as a result of the social distancing intervention, matrix B is replaced by the matrix B[k] in Eq (6); moreover, in order to take into account the effect of vaccination, we assume Δr[k] is subtracted at each day k from the fraction of susceptible individuals, and is simultaneously added to the removed ones, without influencing the fraction of infectious individuals. We reiterate that the effect of the interventions at day k is mediated by the sampling time Δt; in other words, the discrete-time SIR model becomes or, in a compact form

Optimization formulation
The above SIR model with explicit intervention terms is the natural cornerstone for the planning of such interventions.
In particular, we assume a finite-time horizon of k max days and we consider a scenario where at the 0-th day the epidemics is described by given initial conditions s[0], i[0], r[0 The aim of the proposed formulation is to plan the different interventions to be put in place to guarantee the reach of the end of the epidemics on day k max , i.e., we want to enforce dynamical constraints that represent the evolution of the proposed variation of the SIR model (Eq (9), with B[k] and Δr i [k] defined as in Eqs (6) and (7), respectively), together with the requirement that the SIR model reaches the herd immunity surface. The latter requirement is equivalent to enforcing a constraint in the form ensuring that an end-of-epidemic state is reached; at the same time, we want to guarantee that

PLOS ONE
Dynamical intervention planning against COVID-19-like epidemics new infections die out. This requirement, as discussed above, is equivalent to enforcing a constraint in the form of Eq (2). Notably, B[k] is time varying; however, when the final planning instant k max is reached, it is reasonable to assume that non-pharmaceutical interventions are discontinued and E[k max ] = 0 n × n ; thus, the conditions for avoiding epidemic overburst is Let us now discuss the choice variables of the proposed model; specifically, the model aims at identifying the units X[k]�0 n × n and Y[k] � 0 n × n of first and second doses of vaccine to be injected on the k-th day and the intensity of the social distancing measures E[k] 2 [0, e] n × n on the k-th day, for all k 2 {1, . . ., k max }. Notice that, as discussed above, the latter variables must satisfy Let us now focus on aspects related to vaccination. In order to plan for such intervention, we consider a situation where vaccines become available in batches. Specifically, we assume there are specific days k 1 ; . . . ; k w max in which batches of vaccines are received, and we use q½k� 2 R m to denote the vector collecting the total units of vaccines received as of day k for each type of vaccine.
In order to guarantee that the vaccination plan is sound, we need to impose that the cumulative units of vaccine that are injected as of day k do not trespass the received ones, for each type, i.e., Notice that, in order to guarantee that second doses do not trespass the first ones, we consider the constraint in Eq (8); moreover, to guarantee that the overall amount of doses does not exceed the population in each class, we consider a constraint in the form Finally, let us assume that a maximum overall number l h of daily inpatient beds, l sh of which being intensive care inpatient beds, are available. In this view, in order to enforce that the amount of hospitalizations and intensive care hospitalizations do not overcome the limits, we consider constraints in the form where the vectors σ h and σ sh collect the hospitalization and severe hospitalization rates for each class, respectively, and diag[n]i[k] is the vector collecting the population of infectious individuals in each class. Let us now discuss the objective function of the proposed formulation. In particular, we aim to minimize the cumulative intensity of the the social distancing measures over the considered time horizon, i.e., Notice that, within any optimization formulation, minimizing the objective function is secondary to constraint satisfaction. Therefore, within the proposed formulation, reaching of the herd immunity and avoiding the collapse of the healthcare system represent a priority with respect to the minimization of the overall intensity of the social intervention. In other words, solutions that have small objective function value but violate the constraints will be deemed unfeasible and will be discarded by any solver. Overall, the proposed formulation consists of the following Mixed Integer Nonlinear Programming (MINLP) problem.
Rs½k max � � 1 n ; In other words, constraints (I)-(III) model the requirement that the fraction of susceptible, infectious and removed individuals evolve according to the proposed SIR model accounting for the interventions in terms of social distancing and vaccination. Constraint (IV) accounts for reaching an end-of-epidemic state, while Constraint (V) guarantees that new infections die out. Constraint (VI) and (VII) guarantee that the amount of used doses of vaccine do not overcome the available ones or the overall population, respectively. Constraint (VIII) enforces that the second doses (if required) are injected after the adequate time window. Constraint (IX) prescribes that E(k) is symmetric, thus implying that the social distancing effort reducing the influence of the i-th class on the j-th one has a specular effect on the influence of the j-th class on the i-th one. Constraints (X) and (XI) guarantee that the regular and intensive care hospitalizations do not overcome the maximum limit. Finally, constraints (XII)-(XV) guarantee the well-posedness of the variables considered in the formulation.
Approximation strategy. Notice that the proposed formulation amounts to a Mixed Integer Nonlinear Programming problem. In particular, we observe that the model requires a nontrivial amount of variables and constraints (i.e., O(max{n/Δt, n 2 , nm}) for each day of planning. Moreover, we observe that the problem is nonconvex (i.e., considering the nonlinear equality constraints corresponding to the SIR model as two inequality constraints, there is no way both are convex). However, since the units of vaccines involved in the planning are expected to be large, it makes sense to attempt to reduce complexity by considering a continuous relaxation, i.e., dropping integrity constraints. However, also in the case of a convex objective function and a continuous relaxation, the problem has high chances to be NP-Hard (e.g., see [46,47]), thus calling for approximated solutions.
In this paper, our strategy to calculate an approximated solution is to resort to an approximated solver. In fact, we observe that s

[�], i[�], r[�], Δr[�] are actually functions of the variables E[k], X[k], Y[k], even though it is nontrivial to express this dependency in a closed form. Therefore, our strategy is to consider only the variables E[k], X[k], Y[k]
and to express the constraints and the objective function in an algorithmic way, resorting to an approximated solver. Specifically, we use the MIDACO optimization software which implements an extension of the evolutionary Ant Colony Optimization meta-heuristic [48] and which has been developed especially for highly non-linear real-world applications. See [49,50] for a focus of the performance of MIDACO software with respect to the state of the art.
Note that the suggested strategy is independent of a particular solver, but the non-convex nature of the optimization problem suggests an evolutionary approach, like genetic algorithms [51]. Furthermore, the dimensionality of the resulting MINLP in the next case study is very large-scale, consisting of 76650 decision variables and 18295 constraints, and therefore requires a solver that can handle such dimensionality. Finally, we point out that, since in our implementation we chose to evaluate the variables s

[�], i[�], r[�], Δr[�] as a function of the variables E[k], X[k], Y[k]
, a positive consequence is that the step size Δt used to discretize the SIR model has no effect on the overall number of choice variables, which is one of the major sources of complexity for the solution of MINLP formulations (e.g., see [52]).

Computational setting
The optimization with MIDACO was conducted on an Intel 1 Xeon 1 CPU E7 2860 @ 2.27GHz. The CPU runtime for the optimization was fixed to five days. All MIDACO parameters were used by their default values, that means that a feasiblity accuracy of 0.001 was used for all individual constraints listed in Eq (16).

Results
In this section, we test the effectiveness of the proposed formulation by considering a case study with realistic parameters consistent with the current COVID-19 pandemics and relative vaccines. Specifically, we focus on Italy and we identify the optimal vaccination policies over a one-year time horizon, considering 15 age classes (see Table 1) and three types of vaccines. Specifically, Table 2 reports the information regarding the efficacy η j , the delay required to appreciate the effect of the first dose t I j , the delay between doses χ j , the efficacy of the second dose Δη j and the delay required to appreciate the effect of the second dose t II j . The fictional vaccines considered mimic real vaccines, and the parameters are estimates based on data in [53][54][55]. Notably, we assume that the three considered types of vaccine are available only in batches, at specific time instants and in limited amount for each batch, as summarized in Table 3. For simplicity, it is assumed that at regime batches reach between six million and nine million of doses per trimester; such figures are consistent with what has been planned and deployed in Italy [56].
Moreover, we consider a scenario where only a small fraction (i.e., 0.01%) of the age class in thee range 35 − 39 years is initially infected and we assume the maximum daily inpatient beds are l h = 1000, while the maximum daily intensive care inpatient beds are l sh = 100.
Notice that, for the sake of simplicity, we allow vaccination for all age groups, even though Italian regulation does not yet allow COVID-19 vaccination under the age of 5.

Parameter tuning
In order to tune our formulation, we consider the country contact matrix K (see Fig 1), as estimated in [57] for Italy. In particular, only physical contacts have been considered. Notice that the element K ij of a contact matrix from [57] can be considered proportional to the probability that an individual in the i-th age class meets an individual in the j-th; thus, B = Λ � K where

PLOS ONE
Λ ij is the probability that a contact between i and j results in an infection. In this case study, we will use a constant Λ ij = λ; analogously, we will use a constant γ.
To fix the parameters, we will consider a basic reproduction number (i.e., the potential number of new infected generated by one case) R 0 = 3-a value that has been estimated for COVID19 in France [58]. Since for an heterogeneous compartmental model of the form of Eq (1) the role of the basic reproduction number is played by k R k [22,43], we can rescale λ to obtain a basic reproduction number equivalent to the observed one:

PLOS ONE
i.e., The other parameters required to tune the proposed model are the rates of hospitalization, of severe hospitalization, and of death; such parameters can be found in the ECDC ninth risk assessment update for COVID-19 in the EU/EEA and the UK [59] and are reported in Table 1.

Experimental results
Let us now discuss the experimental results from the computational point of view. Fig 2 shows the results of MIDACO in terms of objective function value and overall violation of the constraints, plotted against the number of candidate solutions evaluated by MIDACO. As shown by the figure, we observe that a feasible solution is obtained in about 6 × 10 6 evaluations. As for the objective function, we observe that while the solution is infeasible there is a relevant reduction over time; in particular, we reach a steady solution after about 9 × 10 6 evaluations. Overall, these results suggest the reach of a local minimum.
Having discussed the computational aspects, let us now focus on the structure of the found solution.  Fig 6, where the same data is aggregated and smoothed to improve readability) shows that, in the early stages of the planning, due to the scarcity of vaccines, there is a preference for vaccinating individuals in the age range 20-69 years over young and elderly people; notably, such an age range receives a more or less steady amount of vaccines over time. This stems from the fact that, as discussed above, in the proposed formulation the objective of minimizing the intensity of the social distancing is secondary to constraint satisfaction, i.e., less restrictive social distancing measures can be considered only if they allow

PLOS ONE
the reach the herd immunity and prevent the collapse of the healthcare system. In other words, candidate solutions where strict social distancing was released earlier than the identified solution have been discarded by the solver due to some violation of the constraints. Fig 7 shows the evolution of the proposed SIR model accounting for the effect of social distancing and vaccines, as a result of the interventions planned within the found solution. It can be noted that, due to the strict social distancing measures, only a small fraction of the population becomes infected, with a noticeable peak for the age class in the range 10 − 14 years in correspondence to the softening of the lockdown measures (i.e., around day k = 320). Notice that the fraction of susceptible individuals is slowly eroded due to vaccination, while the fraction of removed has a consequent slow growth due to the resulting immunization. Then, around the end of the planning horizon the lockdown is significantly lifted for the age class in the range 10-14 years (from which the peak in the infected fraction of this age class).
Finally, Fig 9 shows the results of the planning in terms of deaths, hospitalization and intensive care occupancy. It can be noted that the solution found corresponds to a situation where the capacity in terms of regular and intensive care beds is not reached, thus avoiding the collapse of the healthcare system.

Conclusion
In this paper we develop a fine-grained model to support the plan for intervention in order to contrast newborn epidemics. The proposed approach is particularly suitable for infections like the ongoing COVID-19 epidemic, characterized by high infection rate and able to pose the healthcare system and society under stress in terms of intensive care occupancy, deaths or economic consequences. Moreover, the approach allows to plan interventions that blend largescale non-pharmaceutical interventions along with the pharmaceutical ones. Specifically, we build up the planning on two two main types of intervention, namely, non-pharmaceutical (essentially social distancing measures) and vaccination. In order to model the effect of such interventions, we develop a variation of the SIR epidemics model with heterogeneous population; specifically, we assume social distancing intensity to reflect into a reduction of the infection rates, while vaccination to have a partial and delayed immunization effect. Notably, we consider several population classes, several vaccines with different efficacy and with partial and delayed effect, the possibility of a second dose, the availability of vaccines in batches, the need of reaching the herd immunity and the requirement to avoid congestions in the healthcare system. Interestingly, besides representing a detailed, day-to-day planning, the proposed approach also provides useful insights from the clinical and policy-making point of view. In fact, the plan identified by the proposed methodology suggests that, initially, the scarcity of vaccines should be faced by enforcing a strict social distancing, and that vaccination priority should be given to the elderly and "middle-age" population over the younger one. The proposed model exhibits a nontrivial degree of complexity, and the identification of efficient approximated ways to solve it represents a challenging task. Yet, the proposed model represents a remarkably descriptive framework, and future work will be mainly devoted to incorporate other important perspectives for policy and decision makers, such as geographical [61], economical [62] and logistic aspects [63], social equity in the vaccine distribution [64], and skepticism of the population towards vaccines [65].
A last envisaged research perspective is related to the time duration of the planning. In fact, due to the change of the overall epidemiological or pharmaceutical situation, a yearly time horizon could be deemed excessively long; yet, the reach of the herd immunity requires a sufficiently wide time frame. To this end, a viable future work direction is to adopt a "receding horizon" perspective [41,66], where the model is updated after a given period of time (e.g., one month or three months) and a new planning is executed starting from the epidemiological situation at that time. In particular, as new evidence regarding the prevalence of new strains is gathered, the model could be updated (e.g., changing the basic reproduction number [67,68], adding new compartments such as the fraction of asymptomatic individuals [69,70], considering re-infections [71], etc.). Moreover, as new vaccines are developed (or discontinued, as in the case of the Vaxzevria vaccine in Italy [72]) and their effectiveness is better assessed with respect to the circulating variants of the virus, the model can be updated accordingly (e.g., requirement of a booster dose, change in effectiveness, change in the time between doses, etc.). In other words, while the proposed methodology could be considered an open loop approach, we foresee its extension to a closed loop approach. This, of course, raises interesting research questions about the trade-off between the frequency of the update and the computational demands that will be addressed in future work.