Skip to main content
Advertisement
Browse Subject Areas
?

Click through the PLOS taxonomy to find articles in your field.

For more information about PLOS Subject Areas, click here.

  • Loading metrics

Dynamical intervention planning against COVID-19-like epidemics

  • Gabriele Oliva ,

    Roles Conceptualization, Investigation, Methodology, Visualization, Writing – original draft

    g.oliva@unicampus.it (GO); antonio.scala@cnr.it (AS)

    Affiliation Unit of Automatic Control, Department of Engineering, Università Campus Bio-Medico di Roma, Rome, Italy

  • Martin Schlueter,

    Roles Data curation, Software, Validation, Writing – review & editing

    Affiliation Information Initiative Center, Hokkaido University, Sapporo, Japan

  • Masaharu Munetomo,

    Roles Supervision, Validation, Writing – review & editing

    Affiliation Information Initiative Center, Hokkaido University, Sapporo, Japan

  • Antonio Scala

    Roles Conceptualization, Data curation, Investigation, Methodology, Validation, Writing – original draft, Writing – review & editing

    g.oliva@unicampus.it (GO); antonio.scala@cnr.it (AS)

    Affiliations CNR-ISC, Applico Lab, Roma, Italy, Centro Ricerche Enrico Fermi, Roma, Italy, Big Data in Health Society, Roma, Italy, Global Health Security Agenda - GHSA, Roma, Italy

Abstract

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 [25]. 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 [610]; 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 [1820].

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 [2931] and optimization [29, 3235] 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.

Materials and methods

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 Aij. We represent by 0n and 1n vectors with n components, all equal to zero and to one, respectively. Moreover, we use 0n×m, 1n×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 we denote by diag[x] the n × n diagonal matrix having diag[x]ii = xi, for all i ∈ {1, …, n}. On the same footings, given a matrix we denote by diag[A] the n dimensional vector having diag[A]i = Aii, for all i ∈ {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 = AB is such that Cij = Aij Bij; analogously, the Hadamard product between c = ab between two vectors a,b is such that ci = ai bi. Remember that Hadamard product is commutative; moreover, notice that ab = diag[a]b = diag[b]a.

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 . Moreover, let us indicate with s[t], i[t], r[t] the fraction of susceptible, infectious and removed individuals in the -th class at time t and with s[t], i[t], the stack of such variables for all classes. In the following, we assume that s[0], i[0], r[0] ∈ [0, 1]n and s[0] + i[0] + r[0] = 1n. The SIR equations for such an heterogeneous population are given by (1) where B is the n × n transmission matrix, Bij being the rate at which a susceptible individual of class i meets an infectious individual of class j and becomes infected, while the vector collects the rates γi at which infectious individuals in the i-th class are removed from the infection cycle. Notably, the above choice for the initial conditions guarantees that s[t], i[t], r[t] ∈ [0, 1]n and s[t] + i[t] + r[t] = 1n for all time instants t.

End of the epidemics.

Within the SIR model, an epidemic ends when i[t] = 0n, i.e., when such states are also called end-of-epidemic states. At this point, let us define n = [N1, …, Nn]T and let us consider the overall amount of infectious individuals I[t] = iT[t]n; we have that

In this view, since i[t] ≥ 0, the total number of infected individuals I[t] is guaranteed to decrease with time irrespective of the particular value of i[t] if and only if BT diag[n]s[t] − diag[n]γ < 0n, i.e., if and only if [35, 44] (2) where is linked to the next generation matrix [43] that characterizes the stability of an end-of-epidemic state respect to infections and allows to calculate the basic reproduction number R0 indicating the theoretical rate of new infections that an infectious individual could generate.

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 a[t] = i[t] + e[t] 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 with constant step size Δt would be in the form (3)

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 thus obtaining

In other words, we approximate the continuous-time SIR model in Eq (1) via the following discrete-time equations (4) 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 = kt, while we point out that the day corresponding to the time instant hΔt is given by ⌊hΔt⌋.

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 Bij. 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] ∈ [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., Eℓj[k] = 0 means no measure is implemented, while Eℓj[k] = e means the maximum effort is spent in avoiding contacts between the -th and j-th classes. Notice that, by definition, we have that Eℓj[k] = Ejℓ[k]. The above social distancing intensity coefficients account for the different strategies put in place. For instance, if the classes represent geographical regions, then a large value of Eℓj[k] implies large limitation of moving from the -th to the j-th one, while intermediate values of Eℓj[k] could be used to model a scenario where mobility is permitted for work and health circumstances. Conversely, if we consider age classes, then a situation where Eℓj[k] is large for all j could be used to model age-targeted lockdowns (e.g., for the elderly people).

As a result of the choice of Eℓj[k], we consider time-varying terms Bℓj[k] with the following structure (5) or, in a compact form (6) where the n × n matrix E[k] collects the entries Eℓj[k]. In other words, Bℓj[k] corresponds to the nominal Bℓj when no intervention is implemented and reaches zero in the case of a complete lockdown.

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 ∈ [0, 1] after a single dose, while after the second dose the efficacy rises to ηj + Δηj, with Δηj ∈ [0, 1 − ηj]. Notice that the second dose is not required for all types of vaccines; we model this aspect by resorting to a coefficient

Moreover, for each type j of vaccine we assume a time window of days is required for the vaccine to take effect after the first dose, while we use to denote the time window between the first and the second dose (if required) and 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 and a complete effect on day , 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 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 days before with the first dose of the j-th vaccine, weighted by its efficacy ηj, plus the fraction of individuals that were vaccinated 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 (7) where, for the sake of consistency, we assume X[⋅], Y[⋅] are zero when their argument is negative.

Finally, we assume that (8) 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 discrete-time 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 (9) or, in a compact form (10) where z[⋅] = [sT[⋅], iT[⋅], rT[⋅]]T.

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 kmax 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] ∈ [0, 1]n with s[0] + i[0] + r[0] = 1n.

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 kmax, 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 Δri[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 (11) ensuring that an end-of-epidemic state is reached; at the same time, we want to guarantee that 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 kmax is reached, it is reasonable to assume that non-pharmaceutical interventions are discontinued and E[kmax] = 0n × n; thus, the conditions for avoiding epidemic overburst is (12)

Let us now discuss the choice variables of the proposed model; specifically, the model aims at identifying the units X[k]≥0n × n and Y[k] ≥ 0n × 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] ∈ [0, e]n × n on the k-th day, for all k ∈ {1, …, kmax}. Notice that, as discussed above, the latter variables must satisfy (13)

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 in which batches of vaccines are received, and we use 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., (14)

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 (15)

Finally, let us assume that a maximum overall number lh of daily inpatient beds, lsh 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 and 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. (16)

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{nt, n2, 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®Xeon®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 , the delay between doses χj, the efficacy of the second dose Δηj and the delay required to appreciate the effect of the second dose . The fictional vaccines considered mimic real vaccines, and the parameters are estimates based on data in [5355]. 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].

thumbnail
Table 1. Population in the different age classes (Source: [60]).

COVID-19 hospitalization, severe hospitalization and death rates as of April 2021 (source: [59]).

https://doi.org/10.1371/journal.pone.0269830.t001

thumbnail
Table 2. Efficacy of the vaccines considered in the proposed case study.

Vaccine A mimics BNT162b2 (Pfizer & BioNTech); vaccine B mRNA-1273 (Moderna) and vaccine C ChAdOx1 nCoV-2019 (University of Oxford/AstraZeneca). The source for the estimates are: [5355].

https://doi.org/10.1371/journal.pone.0269830.t002

thumbnail
Table 3. Units of Vaccines of each type that are assumed to become available in batches at specific days.

https://doi.org/10.1371/journal.pone.0269830.t003

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 lh = 1000, while the maximum daily intensive care inpatient beds are lsh = 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 Kij 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 Λ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 γ.

thumbnail
Fig 1. Elements of the matrix K of physical contacts among age classes in Italy (source: [57]).

For the sake of readability, the colors of the cells corresponds to ln(Kij).

https://doi.org/10.1371/journal.pone.0269830.g001

To fix the parameters, we will consider a basic reproduction number (i.e., the potential number of new infected generated by one case) R0 = 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 [22, 43], we can rescale λ to obtain a basic reproduction number equivalent to the observed one: i.e., (17)

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 × 106 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 × 106 evaluations. Overall, these results suggest the reach of a local minimum.

thumbnail
Fig 2. Objective function value and overall constraint violation for the solution found using MIDACO, plotted against the number of candidate solutions evaluated.

https://doi.org/10.1371/journal.pone.0269830.g002

Having discussed the computational aspects, let us now focus on the structure of the found solution.

Figs 36 report the structure of the interventions encoded by the found solution. Specifically, according to Fig 3, it can be noted that the the social distancing measures are initially quite intense, and only at the end of the planning horizon there is a partial reduction. Fig 4 shows how vaccine usage is distributed based on the type of vaccine. According to the figure, there is no noticeable difference; this is likely the effect of the scarcity of vaccines in our scenario. Moreover, Fig 5 (as well as 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 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.

thumbnail
Fig 3. Intervention plan corresponding to the found solution in terms of the intensity of social distancing measures plotted against time.

https://doi.org/10.1371/journal.pone.0269830.g003

thumbnail
Fig 4. Intervention plan corresponding to the found solution in terms of the units of the different types administered for each day.

https://doi.org/10.1371/journal.pone.0269830.g004

thumbnail
Fig 5. Intervention plan corresponding to the found solution in terms of the units of vaccines administered to the different age classes for each day.

https://doi.org/10.1371/journal.pone.0269830.g005

thumbnail
Fig 6. Intervention plan corresponding to the found solution in terms of the units of vaccines administered to the different age classes for each day.

The plot aggregates the age classes into the young (0–19 years), middle-age (19–69 years) and elderly (≥70 years) macro-classes. To improve readability, data has been smoothed using a 30-day moving average filter.

https://doi.org/10.1371/journal.pone.0269830.g006

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.

thumbnail
Fig 7. Evolution of the proposed SIR model accounting for the effect of social distancing and vaccines, based on the found solution.

The first, second and third row of plots correspond to the fraction of susceptible, infected and removed individuals, respectively, while the k-th column of plots corresponds to the k-th age class.

https://doi.org/10.1371/journal.pone.0269830.g007

Fig 8 breaks down the social distancing measures by age class; it can be noted that most of the considered time horizon all age classes are strongly restrained in their interaction. 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).

thumbnail
Fig 8. Intensity of lockdown within the found solution for the different age classes and for selected days over the considered time horizon.

The intensity is shown with a blue to yellow scale, where blue represents no social distancing and yellow a complete lockdown.

https://doi.org/10.1371/journal.pone.0269830.g008

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.

thumbnail
Fig 9. Deaths, hospitalizations and intensive care occupancies corresponding to the found solutions, plotted for each day.

https://doi.org/10.1371/journal.pone.0269830.g009

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 large-scale 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.

References

  1. 1. COVID-19 map—John Hopkins Univeristy;. https://coronavirus.jhu.edu/map.html.
  2. 2. Fauci AS, Lane HC, Redfield RR. Covid-19—navigating the uncharted; 2020.
  3. 3. Lau H, Khosrawipour V, Kocbach P, Mikolajczyk A, Schubert J, Bania J, et al. The positive impact of lockdown in Wuhan on containing the COVID-19 outbreak in China. Journal of travel medicine. 2020;. pmid:32181488
  4. 4. Prather KA, Wang CC, Schooley RT. Reducing transmission of SARS-CoV-2. Science. 2020;368(6498):1422–1424. pmid:32461212
  5. 5. Bazant MZ, Bush JW. A guideline to limit indoor airborne transmission of COVID-19. Proceedings of the National Academy of Sciences. 2021;118(17).
  6. 6. Le TT, Andreadakis Z, Kumar A, Román RG, Tollefsen S, Saville M, et al. The COVID-19 vaccine development landscape. Nat Rev Drug Discov. 2020;19:305–306.
  7. 7. Chukwudozie OS, Duru VC, Ndiribe CC, Aborode AT, Oyebanji VO, Emikpe BO. The relevance of bioinformatics applications in the discovery of vaccine candidates and potential drugs for COVID-19 treatment. Bioinformatics and Biology Insights. 2021;15:11779322211002168. pmid:33795932
  8. 8. Cleve M. What the lightning-fast quest for Covid vaccines means for other diseases. Nature. 2021;589:16–8.
  9. 9. Kim JH, Marks F, Clemens JD. Looking beyond COVID-19 vaccine phase 3 trials. Nature medicine. 2021;27(2):205–211. pmid:33469205
  10. 10. Umakanthan S, Chattu VK, Ranade AV, Das D, Basavarajegowda A, Bukelo M. A rapid review of recent advances in diagnosis, treatment and vaccination for COVID-19. AIMS Public Health. 2021;8(1):137. pmid:33575413
  11. 11. Keshavarzi Arshadi A, Webb J, Salem M, Cruz E, Calad-Thomson S, Ghadirian N, et al. Artificial intelligence for COVID-19 drug discovery and vaccine development. Frontiers in Artificial Intelligence. 2020;3:65.
  12. 12. Francis AI, Ghany S, Gilkes T, Umakanthan S. Review of COVID-19 vaccine subtypes, efficacy and geographical distributions. Postgraduate Medical Journal. 2022;98(1159):389–394. pmid:34362856
  13. 13. Beigel JH, Tomashek KM, Dodd LE, Mehta AK, Zingman BS, Kalil AC, et al. Remdesivir for the treatment of Covid-19—preliminary report. New England Journal of Medicine. 2020;383(19):1813–1836. pmid:32445440
  14. 14. Norrie JD. Remdesivir for COVID-19: challenges of underpowered studies. The Lancet. 2020;395(10236):1525–1527. pmid:32423580
  15. 15. Smith T, Killeen GF, Maire N, Ross A, Molineaux L, Tediosi F, et al. Mathematical modeling of the impact of malaria vaccines on the clinical epidemiology and natural history of Plasmodium falciparum malaria: Overview. The American journal of tropical medicine and hygiene. 2006;75(2_suppl):1–10. pmid:16931810
  16. 16. OpenMalaria Software; 2008. https://github.com/SwissTPH/openmalaria.
  17. 17. Van der Ploeg CP, Van Vliet C, De Vlas SJ, Ndinya-Achola JO, Fransen L, Van Oortmarssen GJ, et al. STDSIM: a microsimulation model for decision support in STD control. Interfaces. 1998;28(3):84–100.
  18. 18. van Vliet C, Meester EI, Korenromp EL, Singer B, Bakker R, Habbema JDF. Focusing strategies of condom use against HIV in different behavioural settings: an evaluation based on a simulation model. Bulletin of the World Health Organization. 2001;79:442–454. pmid:11417040
  19. 19. Stuckey EM, Stevenson JC, Cooke MK, Owaga C, Marube E, Oando G, et al. Simulation of malaria epidemiology and control in the highlands of western Kenya. Malaria journal. 2012;11(1):1–14. pmid:23107070
  20. 20. Steen R, Hontelez JA, Veraart A, White RG, de Vlas SJ. Looking upstream to prevent HIV transmission: can interventions with sex workers alter the course of HIV epidemics in Africa as they did in Asia? Aids. 2014;28(6):891–899. pmid:24401648
  21. 21. Kermack WO, McKendrick AG. A contribution to the mathematical theory of epidemics. Proceedings of the royal society of london Series A, Containing papers of a mathematical and physical character. 1927;115(772):700–721.
  22. 22. Diekmann O, Heesterbeek H, Britton T. Mathematical tools for understanding infectious disease dynamics. In: Mathematical Tools for Understanding Infectious Disease Dynamics. Princeton University Press; 2012.
  23. 23. Chen Q. Application of SIR model in forecasting and analyzing for SARS. Beijing da xue xue bao Yi xue ban = Journal of Peking University Health sciences. 2003;35:75–80. pmid:12914225
  24. 24. Laguzet L, Turinici G. Individual vaccination as Nash equilibrium in a SIR model with application to the 2009–2010 influenza A (H1N1) epidemic in France. Bulletin of Mathematical Biology. 2015;77(10):1955–1984.
  25. 25. Bjørnstad ON, Finkenstädt BF, Grenfell BT. Dynamics of measles epidemics: estimating scaling of transmission rates using a time series SIR model. Ecological monographs. 2002;72(2):169–184.
  26. 26. Kühnert D, Stadler T, Vaughan TG, Drummond AJ. Simultaneous reconstruction of evolutionary history and epidemiological dynamics from viral sequences with the birth–death SIR model. Journal of the Royal Society Interface. 2014;11(94):20131106. pmid:24573331
  27. 27. Kretzschmar M, Wallinga J. Mathematical models in infectious disease epidemiology. In: Modern infectious disease epidemiology. Springer; 2009. p. 209–221.
  28. 28. Scala A, Flori A, Spelta A, Brugnoli E, Cinelli M, Quattrociocchi W, et al. Time, space and social interactions: exit mechanisms for the Covid-19 epidemics. Scientific reports. 2020;10(1):1–12.
  29. 29. Behncke H. Optimal control of deterministic epidemics. Optimal control applications and methods. 2000;21(6):269–285.
  30. 30. Nowzari C, Preciado VM, Pappas GJ. Analysis and control of epidemics: A survey of spreading processes on complex networks. IEEE Control Systems Magazine. 2016;36(1):26–46.
  31. 31. Liu F, Buss M. Optimal control for heterogeneous node-based information epidemics over social networks. IEEE Transactions on Control of Network Systems. 2020;7(3):1115–1126.
  32. 32. Duijzer LE, van Jaarsveld WL, Wallinga J, Dekker R. Dose-optimal vaccine allocation over multiple populations. Production and Operations Management. 2018;27(1):143–159. pmid:32327917
  33. 33. Bubar KM, Reinholt K, Kissler SM, Lipsitch M, Cobey S, Grad YH, et al. Model-informed COVID-19 vaccine prioritization strategies by age and serostatus. Science. 2021;371(6532):916–921. pmid:33479118
  34. 34. Enayati S, Özaltın OY. Optimal influenza vaccine distribution with equity. European Journal of Operational Research. 2020;283(2):714–725.
  35. 35. Scala A, Oliva G. Optimal vaccination based on simple, yet effective, linear constraints and vaccines with different effectiveness. Socio-Economic Planning Sciences. Submitted;.
  36. 36. Buckner JH, Chowell G, Springborn MR. Dynamic prioritization of COVID-19 vaccines when social distancing is limited for essential workers. Proceedings of the National Academy of Sciences. 2021;118(16). pmid:33811185
  37. 37. Matrajt L, Eaton J, Leung T, Brown ER. Vaccine optimization for COVID-19: Who to vaccinate first? Science Advances. 2021;7(6):eabf1374. pmid:33536223
  38. 38. Sinha P, Kumar S, Chandra C. Strategies for ensuring required service level for COVID-19 herd immunity in Indian vaccine supply chain. European journal of operational research. 2021;.
  39. 39. Nagurney A. Supply chain game theory network modeling under labor constraints: Applications to the Covid-19 pandemic. European Journal of Operational Research. 2021;293(3):880–891.
  40. 40. Moore S, Hill EM, Tildesley MJ, Dyson L, Keeling MJ. Vaccination and non-pharmaceutical interventions for COVID-19: a mathematical modelling study. The Lancet Infectious Diseases. 2021;21(6):793–802. pmid:33743847
  41. 41. Parino F, Zino L, Calafiore GC, Rizzo A. A model predictive control approach to optimally devise a two-dose vaccination rollout: A case study on COVID-19 in Italy. International journal of robust and nonlinear control. 2021;. pmid:34908815
  42. 42. Burer S, Letchford AN. Non-convex mixed-integer nonlinear programming: A survey. Surveys in Operations Research and Management Science. 2012;17(2):97–106.
  43. 43. Brauer F, Castillo-Chavez C, Castillo-Chavez C. Mathematical models in population biology and epidemiology. vol. 2. Springer; 2012.
  44. 44. Oliva G, Bonfigli S, Scala A. The geometry of herd immunity. In preparation;.
  45. 45. Hall VJ, Foulkes S, Charlett A, Atti A, Monk EJ, Simmons R, et al. SARS-CoV-2 infection rates of antibody-positive compared with antibody-negative health-care workers in England: a large, multicentre, prospective cohort study (SIREN). The Lancet. 2021;397(10283):1459–1469.
  46. 46. Sherali HD, Adams WP. A reformulation-linearization technique for solving discrete and continuous nonconvex problems. vol. 31. Springer Science & Business Media; 2013.
  47. 47. Tawarmalani M, Sahinidis NV. Convexification and global optimization in continuous and mixed-integer nonlinear programming: theory, algorithms, software, and applications. vol. 65. Springer Science & Business Media; 2013.
  48. 48. Dorigo M, Stützle T. Ant colony optimization: overview and recent advances. Handbook of metaheuristics. 2019; p. 311–351.
  49. 49. Schlüter M, Gerdts M, Rückmann JJ. A numerical study of MIDACO on 100 MINLP benchmarks. Optimization. 2012;61(7):873–900.
  50. 50. Schlueter M, Erb SO, Gerdts M, Kemble S, Rückmann JJ. MIDACO on MINLP space applications. Advances in Space Research. 2013;51(7):1116–1131.
  51. 51. Deb K, Pratap A, Agarwal S, Meyarivan T. A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE transactions on evolutionary computation. 2002;6(2):182–197.
  52. 52. Kronqvist J, Bernal DE, Lundell A, Grossmann IE. A review and comparison of solvers for convex MINLP. Optimization and Engineering. 2019;20(2):397–455.
  53. 53. Polack FP, Thomas SJ, Kitchin N, Absalon J, Gurtman A, Lockhart S, et al. Safety and efficacy of the BNT162b2 mRNA Covid-19 vaccine. New England Journal of Medicine. 2020;. pmid:33301246
  54. 54. Baden LR, El Sahly HM, Essink B, Kotloff K, Frey S, Novak R, et al. Efficacy and safety of the mRNA-1273 SARS-CoV-2 vaccine. New England journal of medicine. 2021;384(5):403–416. pmid:33378609
  55. 55. Knoll MD, Wonodi C. Oxford–AstraZeneca COVID-19 vaccine efficacy. The Lancet. 2021;397(10269):72–74.
  56. 56. OpenData on Italian Vaccines;. https://github.com/italia/covid19-opendata-vaccini.
  57. 57. Mossong J, Hens N, Jit M, Beutels P, Auranen K, Mikolajczyk R, et al. Supporting materials for Social contacts and mixing patterns relevant to the spread of infectious diseases. PLoS medicine. 2008;5(3):e74.
  58. 58. Salje H, Tran Kiem C, Lefrancq N, Courtejoie N, Bosetti P, Paireau J, et al. Estimating the burden of SARS-CoV-2 in France. Science. 2020;369(6500):208–211.
  59. 59. Rapid Risk Assessment: Coronavirus disease 2019 (COVID-19) in the EU/EEA and the UK– ninth update;. https://bit.ly/39Kz9FO
  60. 60. US Census Bureau;. https://www.census.gov/.
  61. 61. Galeazzi A, Cinelli M, Bonaccorsi G, Pierri F, Schmidt AL, Scala A, et al. Human mobility in response to COVID-19 in France, Italy and UK. Scientific reports. 2021;11(1):1–10. pmid:34162933
  62. 62. Bonaccorsi G, Pierri F, Cinelli M, Flori A, Galeazzi A, Porcelli F, et al. Economic and social consequences of human mobility restrictions under COVID-19. Proceedings of the National Academy of Sciences. 2020;117(27):15530–15535.
  63. 63. Organization WH, et al. COVID-19 vaccination: supply and logistics guidance: interim guidance, 12 February 2021. World Health Organization; 2021.
  64. 64. Wu JH, John SD, Adashi EY. Allocating vaccines in a pandemic: the ethical dimension. The American Journal of Medicine. 2020;133(11):1241–1242. pmid:32653419
  65. 65. Briand SC, Cinelli M, Nguyen T, Lewis R, Prybylski D, Valensise CM, et al. Infodemics: A new challenge for public health. Cell. 2021;184(25):6010–6014. pmid:34890548
  66. 66. Agachi PS, Cristea MV, Csavdari AA, Szilagyi B. 2. Model predictive control. In: Advanced Process Engineering Control. De Gruyter; 2016. p. 32–74.
  67. 67. Ito K, Piantham C, Nishiura H. Relative instantaneous reproduction number of Omicron SARS-CoV-2 variant with respect to the Delta variant in Denmark. Journal of medical virology. 2022;94(5):2265–2268. pmid:34967453
  68. 68. Ito K, Piantham C, Nishiura H. Predicted dominance of variant Delta of SARS-CoV-2 before Tokyo olympic games, Japan, July 2021. Eurosurveillance. 2021;26(27):2100570.
  69. 69. Di Giamberardino P, Iacoviello D, Papa F, Sinisgalli C. Dynamical evolution of COVID-19 in Italy with an evaluation of the size of the asymptomatic infective population. IEEE Journal of Biomedical and Health Informatics. 2020;25(4):1326–1332.
  70. 70. Kheifetz Y, Kirsten H, Scholz M. On the parametrization of epidemiologic models–lessons from modelling COVID-19 epidemic. arXiv preprint arXiv:210911916. 2021;.
  71. 71. Wilhelm A, Widera M, Grikscheit K, Toptan T, Schenk B, Pallas C, et al. Reduced neutralization of SARS-CoV-2 omicron variant by vaccine sera and monoclonal antibodies. MedRxiv. 2021;.
  72. 72. European Centre for Disease Prevention EC, Control. Overview of EU/EEA country recommendations on COVID-19 vaccination with Vaxzevria, and a scoping review of evidence to guide decision-making. 2021.