Epidemic control by social distancing and vaccination: Optimal strategies and remarks on the COVID-19 Italian response policy

1.   Introduction

The dramatic worldwide impact of the COVID-19 pandemic and the multiple failures in the related responses ^[1] have highlighted a number of challenges for future preparedness activities. Among these challenges, three have emerged as the most critical ones. The first one is represented by the trade-off between the direct and the indirect costs of the epidemic, that is, the tension between the direct effects of an outbreak (in terms of deaths and overwhelming of public health resources) and its indirect impacts –at the broader societal level– resulting from the undertaken control measures. This leads to the second major challenge, namely which are the appropriate measures for balancing these costs' trade-off. The COVID-19 debate has highlighted a number of different control options including mitigation, suppression, aggressive containment, and elimination ^[2,3,4,5]. The last challenge concerns the development of fully agreed principles for unambiguously ranking different control options in preparedness protocols for measuring the aforementioned costs' trade-off, which appears to be a challenge per se ^[6].

Optimal control potentially appears to be a critical tool for the development of the aforementioned principles. The COVID-19 pandemic has seen a true explosion of contributions applying optimal control methodologies to epidemic phenomena and their various facets, including the optimal implementation of non-pharmaceutical interventions (NPI), such as social distancing, testing tracing and isolation of infective cases as well as pharmaceutical ones, or vaccine prioritization (see e.g., ^{[6,7,8,9,10,11,12,13,14,15,16]} and references therein).

This work aims to contribute to this area by investigating the optimal control of a threatening epidemic of an age-structured infection. Consistently with the COVID-19 experience, we seek the optimal switch conditions between a control regime only relying on non-pharmaceutical interventions (NPI, such as social distancing) in a first epoch of the epidemic, to a new control regime in a subsequent epoch, based on an appropriate combination of social distancing and vaccination, thanks to the arrival of an effective vaccine allowing a gradual relaxation of social distancing. With this aim, we combine and extend a number of previous works. The epidemiological framework considers a parametrized SIRS model for a communicable disease, structured by age of infection, and used for analysing multi-phasic epidemics ^[17,18]. The model relies on the hypothesis of a high prioritization of health protection from the epidemic attack, that we term the low attack rate (LAR) hypothesis. The LAR hypothesis assumes that, regardless of the epidemic phase, an appropriate degree of social distancing is always mantained to prevent the susceptible fraction from being substantially depleted by epidemic incidence. This hypothesis was proven possible in all countries worldwide that opted for elimination, e.g. China, Denmark, New Zealand, Southern Korea, and Australia ^[3,4,5]. The proposed model was plugged into a theoretical, finite-horizon, optimal control setting with NPIs and vaccination ^[19], where a number of formal results on existence, uniqueness, and basic properties of the optimal action were determined.

In particular, in this article, we develop the cited works in order to implement the optimal control solution. To do so, we first provide a detailed formulation of the various components of costs raised by the epidemic (direct, indirect, vaccination). Next, we propose a specific numerical algorithm for solving the optimal control problem ^[19] via a discretization of the forward state equations and the corresponding dual (backward) problem. The proposed algorithm is used to analyse the patterns of the optimal strategy. Amongst other things, we discuss the dependency of the optimal solution on the degree of prioritization on the direct costs of the epidemic and we compare the optimal solution with actual data from the COVID-19 experience in Italy. To do so, we parametrize the components of the cost function by combining a range of literature data and studies.

The manuscript is organized as follows: Section 2 reports the general epidemiological model; Section 3 describes the costs components and the optimal control problem; Section 4 reports the general results of the optimization procedure; the parametrization of the problem in the form of differential delay equations (DDE) is reported in Section 5, while the corresponding numerical algorithm is exposed in Section 6; the results are presented in Section 7; the discussion and concluding remarks follow.

2.   The epidemiological model

Following ^[19], we consider a SIRS model of a threatening outbreak of a communicable infection in a constant population. The model yields temporary immunity (both natural and vaccine-related) and it is structured by age-of-infection (sometimes called time-since-infection) along the tradition of the basic epidemic model by Kermack and McKendrick ^[20]. Following this model, we also assume that the epidemic time-scale is short with respect to demographic processes so that we can disregard the vital dynamics of the population.

As has been the case for COVID-19, the outbreak is assumed to be controllable by two main types of interventions, namely

$({\rm{i}})$ non-pharmaceutical interventions, primarily represented by social distancing and mainly affecting the transmission rate over time,

$({\rm{ii}})$ vaccination.

As expected for pandemic attacks, social distancing is assumed to be the only available intervention in a first phase, while a combination of social distancing and immunization is considered after a vaccine becomes available. For the sake of simplicity, we disregard testing/tracing and isolation of infective individuals, which can be incorporated into the removal process. Moreover, we will rely on the LAR hypothesis mentioned in the introduction ^[17,19]. The LAR hypothesis postulates that, after an initial phase of epidemic growth, effective control is undertaken, such that:

$({\rm{i}})$ during the pre-vaccination period (Period 1 from now on), early and effective social distancing measures are enacted, which adequately lower infection transmission. This allows a slow depletion of the susceptible population so that the susceptible fraction remains close to 100%;

$({\rm{ii}})$ after the introduction of vaccination, i.e., during the post-vaccination period (Period 2), the epidemic remains adequately controlled through a combination of social distancing measures and immunization. In this period, vaccination becomes the key responsible for the depletion of the susceptible population.

On these assumptions, when the vaccination campaign starts, the susceptible population is still essentially nearby $100\%$, after the vaccination campaign progresses, the removed population is essentially composed of vaccinated individuals. Under the LAR hypothesis, the model proved effective in describing the first two years of the COVID-19 epidemic in Italy ^[17]. The LAR hypothesis has a number of advantages, both regarding modeling and substantive aspects. First, it can be taken as a crude strategy to prevent epidemic incidence to grow dramatically, which avoids including –at least as a first approximation– external constraints to account for the finite capacity of public health resources (e.g., hospitals). Second, it avoids the need to explicitly include infection-related mortality. By also disregarding vital dynamics, the population can be taken to be constant over time with overall size $N$. From a substantive viewpoint, the LAR hypothesis represents a manner to include a high degree of prioritization by public policy on controlling the direct cost of the epidemic. We recall that, for COVID-19, the LAR hypothesis is far from academic, as it was surely correct for those countries worldwide that opted for elimination or for aggressive suppression rather than mitigation, as has been the case of China, New Zealand, Australia, South-Korea, Japan, etc. ^[3].

2.1.   The state system

Our basic model is described by the following system structured by age of infection ^[17,19]

endowed with the initial condition

In (2.1), ${ } s(t)$ denotes the susceptible fraction, ${ } v(t)$ the per-capita number of daily successful immunizations, and

is the absolute density of infected individuals having an age of infection equal to $x$ at time $t$, i.e., those who acquired the infection $x$ time units before (that is, at time $t-x$). Finally $x_+$ denotes the maximum age of infection. Note that we are not considering demographic dynamics, so that the total population size $N$ is constant overtime.

Model (2.1) includes the LAR hypothesis in the dynamics of the susceptible fraction, which therefore reflects the trends of vaccination and immunity losses only because, with respect to the usual formulation of SIR and SIRS models, the contribution of infection incidence is considered to be negligible. The trend of infective individuals along the age of infection is described by the so called Von Foerster-McKendrick partial differential equation (2.1, $ii$), summarizing the dynamics along characteristic lines, by the related boundary condition (2.1, $iii$), representing the incidence of new infections per unit of time (i.e., infective individuals having age of infection equal to $0$), and by the initial condition (2.2). The model parameters have the following meaning

● $v(t)$ = the per person rate of effective (i.e., successful) immunization per unit of time,

i.e., coarsely speaking, the fraction of newly vaccinated susceptibles that are effectively protected against infection. Function $v(t)$ is identically zero prior to the arrival of the vaccine and has an upper bound $v_{max}$ determined by the capacity of the vaccination system thereafter. In the experience of industrialized countries, which were the first to start mass vaccination campaigns against COVID-19 in view of their market power, the distribution of vaccines over time has shown a complicated range of shapes, reflecting, e.g., logistic difficulties in an early phase, supply shortages, temporary plateauing due to achievement of maximal capacity, vaccine scare phenomena (e.g., due to the VaxZevria) alert, or age-dependent vaccine-hesitant behaviors (see e.g., the discussion in ^[18]).

● $\delta$ = the immunity waning rate,

i.e., the rate at which removed people loose their immunity; for the sake of simplicity, $\delta$ is taken here to describe both natural and vaccine waning immunity. For modeling reasons ^[19], we will assume

● $c(t)$ = the number of adequate contacts per person and per unit of time, where $c(0) = c_0$ is the reference value under normal conditions, i.e., in the absence of the epidemic. We assume:

where $c_m$ is the minimum affordable number of contacts. In the subsequent optimal control application, we will conveniently use as control variable the ratio between the altered level of contacts prevailing under the interventions enacted at time $t$ and its normal level

The previous normalized quantity, that we term the altered contact ratio, is the simplest measure of the contacts' alteration under social distancing interventions.

● ${\beta}_0(x)$ = average intrinsic infectiousness of an infected individual aged $x$,

i.e, the probability that an infected individual aged $x$ infects a susceptible during an adequate contact, possibly related to their viral load. We deliberately did not include time dependencies in the infectiousness due to individual protections measures, e.g., wearing masks, given that we can simply assume that the latter measures scale the adequate number of contacts and are therefore embedded into $c(t)$.

● $\gamma(x)$ = the removal rate of infected individuals aged $x$ at time $t$,

i.e., the overall per-capita rate at which infected individuals aged $x$ are removed by any cause as e.g., recovery or death (during an uncontrolled epidemic) or by screening, tracing, isolation, etc, in the presence of interventions. For simplicity, $\gamma(x)$ is assumed to be independent of time. Note that ${\gamma}(x)$ is related to the probability $\Gamma(x)$ to still be infected after a time $x$ since infection (the so-called survival-to-removal function), defined by

From the previous definitions, we can compute key epidemiological parameters such as the transmission rate of infected individuals aged $x$ at time $t$

the generation time distribution

the basic reproduction number

and the effective reproduction number

In particular, we can express the incidence of infection using ${\rho}(t)$

where

is the incidence in the absence of distancing measures and in a totally susceptible population.

Problem (2.1) provides the state system for the control problem we introduce in next section. However, as it is shown in ^[19], this problem is equivalent to

with

where $Z(t)$ is the variable defined in (2.6) and related to incidence by (2.5). The connection with system (2.1) is provided by the well-know formula for the density of infection prevalence

Note that both quantities $\mathcal{K}$ and $F$ have a sharp epidemiological interpretation. The former represents the number of secondary infections by a typical infective individual with age of infection $x$ when social activity takes place at the altered level $c(\sigma)$ and the susceptible fraction equals $s(\sigma)$. The latter represents the number of secondary infections collectively generated by the initial age distribution of infective individuals.

In (2.7) the couple $({\rho}(t), v(t))$ is assigned and, when necessary, we will denote by $s[v](t)$ and $Z[{\rho}, v](t)$ the solutions to (2.7), showing their dependence on the input.

2.2.   The initial density of infected

Within the proposed framework, a particular attention must be paid to the initial density $Y_0(x)$. This function captures the information embedded in the initial phase of free epidemic growth before the control action is enacted at time $t = 0$. This free epidemic early growth is assumed to take place over a time interval $[-T_0, 0]$ during which no public control measures are undertaken and spontaneous behavior changes are assumed to be negligible. This phase is characterized by

and we have

The latter relation follows from the LAR hypothesis and states that the effective reproduction number ${{\cal R}}_E(t)$ (summarizing epidemic growth per generations during an arbitrary phase of the epidemic) coincides with the the basic reproduction number ${{\cal R}}_0$, summarizing epidemic growth per generations in a wholly susceptible population (and in the absence of intervention measures). Following ^[17], we suppose that this initial infective cohort is distributed according to the stable age distribution (SAD) corresponding to the given epidemiological parameters

where ${\alpha}^*$ is the leading root of the characteristic equation

The previous hypotheses imply that, for $t\in [-T_0, 0]$,

where $I_0$ is the total number of infected individuals $I(t)$ at $t = 0$. In fact we have

where $\alpha^*$ also represents the growth exponent of the infected in the initial phase of the epidemics. Note that the previous hypothesis is indeed quite mild: it just assumes that the phase of free epidemic growth is long enough to allow transients due to the "very initial" epidemic seeds to disappear and the ensuing SAD to emerge.

3.   The cost function

Model (2.1) represents the state system for our optimal control problem. The functions ${\rho}(t)$ and $v(t)$, i.e., the altered contact ratio and the effective immunization rate, will play the role of the control variables designed to minimize an appropriate cost function over a horizon having finite duration $T$. In particular, the date of vaccine arrival ($T_1 < T$) is assumed to be known with certainty. The COVID-19 experience has shown that the latter hypothesis on the timing of vaccine availability is not unreasonable. Indeed, after the full mapping of the viral genome in China, the early predictions on the date of vaccine arrival were revealed to be surprisingly correct, possibly due to the massive inflow of public resources from industrialized countries to the vaccine industry. For early instances of optimal control of COVID-19 under random vaccine arrival time, see the economic contributions in ^[7,8].

As pointed out in the introduction, an appropriate cost function should include both the direct cost of the epidemic, primarily due to deaths and hospitalizations, and its indirect cost, resulting from the societal injuries (i.e., economic, social, relational, and non-epidemic-related health) arising from the measures enacted to control the outbreak. The rationale is that the same pool of social contacts that are responsible for infection transmission are also the medium promoting all structured economic activities, e.g., production and consumption ^[21], as well as any other type of social relation. Consequently, control measures cutting social contacts will have a negative feedback on any type of human activity, thereby decreasing the total societal welfare. Here, we do not enter into the issue of how exactly to measure societal welfare and simply follow the economists' practice to rely on gross domestic product (GDP) as was also done for COVID-19 ^[7]. Therefore, in our evaluation of costs, we will assume that total societal welfare is proportional to GDP. Clearly, this must be considered a first step, as it is widely known that GDP is just a measure of material production and, therefore, that other indicators (e.g., the human development index by the UN) should be considered.

3.1.   Direct cost of the epidemics

As for the direct health costs of the epidemic, namely those resulting from cases of serious disease after infection (primarily hospitalizations and deaths), these arise from the total epidemic incidence over the whole period $[0, T]$,

where $(s(t), Z(t))$ is the solution of the state system (2.7) with inputs ${\rho}(t)$ and $v(t)$, and $c_E$ is the average cost generated by one infected individual per unit of time over period $[0, T]$. Here, we choose the following form for $c_E$

where

● $g \in (0, 1)$ is the fraction of incident cases yielding to serious disease. From (3.2), this means that we do not attribute costs to the fraction $(1-g)$ of non-serious cases.

● $\mu$ and $(1-\mu)$ ($\mu \in (0, 1)$) are the fractions of serious cases who die due to the disease or do not die, respectively. Surviving individuals are assumed to generate a cost including hospitalization costs as well as costs for treating subsequent sequelae (e.g., long-COVID). For simplicity, we do not associate a hospitalization cost to dead individuals.

● $\kappa_H$ is the average monetary cost of an hospitalization for a survivor.

● $\kappa_D$ is the average cost of one individual death expressed in monetary terms for comparing it with other cost items. We apply here a concept of emotional cost i.e., the death of an individual will cancel their entire contribution to welfare at all future times regardless of his age. Therefore, this cost is related to the residual (average) life expectancy $T_L$ of an individual dying due to the disease. In what follows, we will assume that $\kappa_D = T_L \times GDP$, where $GDP$ represents the per-capita income of the community afflicted by the outbreak. The underlying idea is that GDP is taken as a proxy for the welfare flow generated by a living person living during one year, regardless of whether the person is actually working or not. Other choices are clearly possible: for example, economic approaches have often used the value of a statistical life ^[7,8].

3.2.   Indirect cost of the epidemic

To represent indirect costs, which arise from intervention $\rho$ at time $t$, we introduce a loss function $Q({\rho})$ with the following properties

where

● ${\rho} = 1$ corresponds to normal social activity, i.e., when no social distancing measure is in place and the level of social activity stays at its normal level $c_0$;

● ${\rho} = {\rho}_m$ represents the minimum affordable reduction of social activity. This corresponds ideally to a level of social contacts allowing the enactment of those critical activities whose prolonged cancellation would create irreparable damages at the societal level. Clearly, ${\rho}_m$ corresponds to a maximum of indirect costs (and related welfare loss).

Function $Q({\rho})$ is such that

represents the overall societal welfare loss per unit of time when per-capita social contacts are set to the altered level ${\rho}$, ${\rho}\in [{\rho}_m, 1]$. From (3.3), the maximum loss ${{\cal L}}$ occurring at ${\rho}_m$ is clearly approached under harsh lockdown conditions.

A form of function $Q(\cdot)$, which is convenient for our application, is the following

This formulation has a standard quadratic form for $\omega = 0$. As $Q'(1) = \omega$, the parameter $\omega\ge 0$ tunes the reaction to changes of social contacts when departing from the reference value $c_0$. In particular, form (3.4) is normalized, given that $Q(\rho_m) = 1$.

Through the previous function, the overall indirect cost arising from a prescribed social distancing policy ${\rho}(t)\ (t\in [0, T])$ is given by:

In the special case of a constant policy, with ${\rho}(t)$ set to a constant level $\overline{{\rho}}$ for the entire horizon, the overall loss would amount to

The previous definition of the indirect cost of social distancing is a minimalistic one, straightforwardly reflecting the standard epidemiological approach by which any reduction of the social activity from its normal level due to the enacted restrictions, yields –by definition– a variety of societal costs (e.g., ^[15] and references therein). More explicit representations of indirect costs due to lockdowns were adopted already in seminal economic efforts during the first wave of the COVID-19 literature ^[7,8], where such costs were represented through the loss of income production (and consumption) arising from reducing work and production activities. In particular, the quadratic formulation (3.4) represents a simple parametric family of functions that (ⅰ) are not too sensitive to mild social distancing interventions, i.e., for small departures from the normal contact rate ${\rho} = 1$ (corresponding to absence of epidemic alert); and (ⅱ) generate dramatically increasing indirect costs when large alterations to social contact activity are implemented. This reflects the idea that large social distancing policies, breaking social, economic, and relational ties, will cause dramatic costs on the society as a whole.

3.3.   Cost of the vaccination campaign

Finally, the component of total cost due to the vaccination campaign over the period $[0, T]$ is taken as

where $V(t)$ is the absolute vaccination rate and $c_V(v(t))$ is the gross cost associated with the administration of a single vaccination course. This cost goes beyond the material cost of the single vaccine dose paid by the government and includes all the costs associated with the campaign during a pandemic, namely organization and logistic costs, search costs for reaching elusive or hesitant people, increase of costs in the presence of a sustained worldwide demand for vaccines, etc. We assume that this cost depends on the vaccination rate according to the following general form (still expressed in terms of per-capita welfare per unit of time)

where the function $H(\cdot)$ satisfies

In equation (3.8) the multiplicative constant $\eta$ represents the basic cost at "low" vaccination rates, while function $H(\cdot)$ introduces a (weak) growth cost component related to the aforementioned capacity factors that arguably prevent exact proportionality. The simplest form for $H(\cdot)$ is

with $0 < \varepsilon \ll 1$. This is the form we will adopt in our simulations. This nonlinear form reflects the idea that the gross vaccination cost is largely proportional to the number of vaccinations administered during the campaign (obviously including organizational costs, etc, ) but must also include a nonlinear term accounting for the fact that achieving a high vaccination coverage produces additional costs for the system that do not grow linearly, for example due to the difficulties that arise when the system plans to reach marginal or elusive subpopulations of individuals. The latter case can include strongly hesitant people but also less informed groups, such as the foreign-born or low-literacy groups (see, e.g., ^[22]).

3.4.   The overall cost function

Given a control strategy of outbreaks, here represented by a prescribed combination $({\rho}(t), v(t))$ of a social distancing action and an immunization campaign, the social planner will combine the various related cost components ${{\cal C}}_E, {{\cal C}}_w{{\cal C}}_V$ into the overall cost functional. However, as well documented by the COVID-19 experience, the structure to attribute to the overall cost functional is far from straightforward. Indeed, as pinpointed in the Introduction, during the COVID-19 pandemic, countries followed very different control attitudes and strategies ranging from elimination, as followed by China and later on by Australia and New Zealand, to various forms of suppression or mitigation up to the mild control followed by Sweden for the entire first pandemic year. These different options can be seen to reflect a widely different prioritization of the various cost items by the social planner. We incorporate these differences on government-level preference for direct versus indirect costs by attributing a priori weights on costs components. Namely, we take total costs as

where the weights $\chi$ and $(1-\chi)$ represent the relative governmental preference for direct versus indirect cost, and vice versa. A similar representation but without vaccination has been used ^[15]. This formulation merges the direct cost of the infection and vaccination as a unique cost item to be compared with indirect costs. The merging is motivated by traditional approaches to optimal control of vaccine-preventable infectious diseases, where direct infection costs and vaccination costs are the main cost items. Said otherwise, the first component of (3.10) reflects the total direct costs of the policy, i.e., the sum of the direct cost of infection plus the direct cost of the vaccination control on the assumption that the (direct) cost of social distancing (i.e., keeping people isolated at home) is negligible. The second term in (3.10) represents the (weighted) contribution of indirect costs arising from the socio-economic losses caused by the undertaken interventions. In this perspective, (3.10) represents the correct trade-off between the direct and indirect costs of the policy.

In the following, we will adopt the constitutive forms (3.4) (3.9) and will be concerned with minimizing the functional

with

The functional (3.11) is to be minimized in the set of controls ${{\cal U}} \subset L^\infty(0, T)\times L^\infty(0, T)$, defined as

reflecting the partition into the two epidemic periods, Period 1 and Period 2, discussed above.

In its full form, problem (3.11)-(3.12) implies to seek an optimal social distancing schedule during Period 1 ($[0, T_1)$) and an optimal combination of social distancing and vaccination during Period 2 ($[T_1, T)$). We term it an optimal switching problem. The problem encompasses a number of relevant sub-cases, namely

● the case $T_1 = T$, amounting to seek the optimal social distancing action over a finite horizon in the absence of vaccination. This has represented a major task of optimal control research on COVID-19 in an earlier phase but might also represent the situation where the pathogen characteristics rule out the possibility of developing a vaccine in any reasonable time horizon, as it has been the case for the HIV pandemic ^[23].

● the case where the vaccination schedule (during Period 2) is taken as given due to, e.g., external constraints on the vaccine supply side, as happened for COVID-19. This case collapses again into one of optimal social distancing only.

● all relevant sub-cases resulting from attributing different values to the weight of direct (indirect) costs $\chi$.

4.   Determining the optimal strategy

Based on (2.7), in ^[19] we have proved the existence of an optimal control and, in order to state the optimality conditions satisfied by an optimal couple $({\rho}^*(t), v^*(t))$, we have considered the (dual) system coupled to (2.7)

In fact, denoting $p[{\rho}, v](t)$ and $q[{\rho}, v](t)$ the solutions to (4.1), we proved that an optimal couple $({\rho}^*(t), v^*(t))$ must satisfy the system of equations

where $N_{[{\rho}_m, 1]}$ and $N_{[0, v_{max}]}$ are the normal cones of the intervals $[{\rho}_m, 1]$ and $[0, v_{max}]$, respectively.

Now, since the two (multivoque) functions $M_1: {{\bf R}} \to {{\bf R}}$ and $M_2: {{\bf R}} \to {{\bf R}}$ defined as

are strongly maximal monotone, and the solutions $p[{\rho}^*, v^*], q[{\rho}^*, v^*], s[v^*], Z[{\rho}^*, v^*]$ are Lipschitz continuous in $L^\infty(0, T)$, system (4.2) has a (possibly) unique solution for the optimal control actions on social distancing and (effective) vaccine administration. This solution can be obtained through the proximal point algorithm ^[24], namely through the approximating sequence $({\rho}^k, v^k)$ defined by recurrence from any starting point $({\rho}^0, v^0)$:

where ${\lambda}$ is an auxiliary parameter aiming to optimize the algorithm. Indeed, the sequence $({\rho}^k, v^k)$ converges strongly in $L^\infty(0, T)$ to the solution of (4.2). The numerical algorithm we will set up in Section 6, is a discrete version of this basic result.

Note that with the definitions (4.3) we have

and

5.   Reduction to DDE systems via Erlang kernels

In a previous paper ^[17], we adopted a special parametrization for the age-dependent epidemiological functions $\beta_0(x), \gamma(x)$ allowing reduction of problem (2.7) to a (possibly delayed) system of ordinary differential equations. Thus, in order to cast our problem in a DDE framework we assume that both ${\beta}_0(x)$ and $\gamma(x)$ obey translated Erlang-type functions; namely, we take $x_+ = +\infty$ and

${\rm{i}})$ the infectivity kernel ${\beta}_0(x)$ is a (non-proper) translated Erlang density of order $2$ and rate $\varphi$:

Note that it holds

${\rm{ii}})$ the removal rate ${\gamma}(x)$ is taken as

corresponding to the following survival-to-removal function:

These specific functions ${\beta}_0(x)$ and $\gamma(x)$ belong to an extended class of Erlang-type functions allowing reducibility that was considered in ^[17]. The specific minimal choice above was used to parsimoniously describe the pre-vaccination course of the COVD-19 epidemic in Italy ^[17,18]. Here, we will use the same choice to approach the control problem presented in the previous section. Within the context determined by (5.1) and (5.2), the basic reproduction number ${{\cal R}}_0$ reads

where $\theta = {\gamma} + \phi$. Also, we have the characteristic equation (2.10)

yielding the generation time distribution

which, therefore, obeys a $\tau$ translated Erlang density of index (or shape) parameter $2$ and rate parameter $\theta$.

Now, using (5.1) and (5.2) and introducing the auxiliary variable

we can replace (2.7) with the following system

endowed with the initial heredities for $t\in [-\tau, 0]$ (see (2.11))

and the initial condition

where

Indeed (see ^[17]), we can draw (5.6) and (5.7) differentiating $Z(t)$ in (2.7) using the specific forms (5.3) and (5.5). Based on the previous variables, we can compute other significant variables on the time period $[0, T]$ such as the number of infected people $I(t)$ or infective people $I^\#(t)$. Actually, we will use this latter to compare the theoretical result with real data: it will obey to the equation

Also, in order to use the characterization (4.2), we similarly transform the dual system (4.1). Namely, using the auxiliary variable

we can reduce (4.1) to the following (backward) DDE system

endowed with the conditions

Note that the conditions on the interval $[T-\tau, T]$ arise because the kernel $K(x)$ vanishes on the interval $[0, \tau]$; moreover, equations (5.11, $i$) and (5.11, $ii$) follow differentiating the respective integral equations.

The state system and the dual one will be discretized in order to set up a numerical algorithm for calculating an optimal control.

6.   Numerical algorithm

The algorithm we propose provides a discrete approximation of an optimal couple $({\rho}^*, v^*)$ by first discretizing the forward state problem (5.6) and the backward dual problem (5.11) with a first-order numerical method.

Considering a time step $h$, we discretize the time interval $[-\tau, T]$ by defining

Then, the whole interval $[-\tau, T]$ is described for $n = - n_{\tau}, \dots, n_{T}$ and

Concerning the state system (5.6), we consider $S_n$, $Z_n$, $J_n$ as the approximation of $s(t_n)$, $Z(t_n)$, $J(t_n)$, to be computed through the following forward recurrent scheme for $0 \le n \le n_T$

starting with the initial heredities (see (5.7) - (5.9))

and the initial condition

In (6.1) the approximated values of the vaccination rate $v$ and of the contact ratio $\rho$

are assigned and satisfy the constraints

Related to (6.1), we also compute the approximated values of the number of infective individuals $I^\#_n = I^\#(t_n)$ (see (5.10)) by the scheme

Besides, setting $P_n = p(t_n)$, $Q_n = q(t_n)$, $R_n = r(t_n)$, we discretize (5.11) by the backward recurrent scheme, for $0 \le n \le n_T-n_\tau$

initialized by the conditions

while $P_n$ for $(n_T-n_\tau) \le n \le n_T$, is obtained by backward recurrence through

The recurrent schemes (6.1) and (6.3) are the main tool for the algorithm, aiming to compute an approximation of an optimal couple. Note that in (6.1) the vectors $C\equiv (C_{-n_\tau}, \dots, C_{n_T})$ and $V\equiv (V_{0}, \dots, V_{n_T})$ are assigned, while in (6.3) the assigned input is given by the vectors $C\equiv (C_{-n_\tau}, \dots, C_{n_T})$, $S\equiv (S_{-n_\tau}, \dots, S_{n_T})$ and $Z\equiv (Z_{-n_\tau}, \dots, Z_{n_T})$, with $S$ and $Z$ coming from (6.1). Thus, in view of our algorithm, we will denote

the solutions of (6.1) and (6.3), respectively, showing the respective inputs. Then, we propose the following procedure based on the proximal point algorithm (4.4):

${\bf{1}}$ choose $0 < {\lambda} \le 1$ and a starter couple $C^0\equiv (C^0_{-n_\tau}, \dots, C^0_{n_T})$, $V^0\equiv (V^0_{0}, \dots, V^0_{n_T})$ satisfying (6.2). Then, for $k\ge 0$, follow the iterative procedure from step 2

${\bf{2}}$ compute ${{\cal S}}(C^k, V^k)$ via (6.1), providing $S^k$, $Z^k$ and $J^k$

${\bf{3}}$ compute ${{\cal D}}(C^k, S^k, Z^k)$ via (6.3), providing $P^k$, $Q^k$ and $R^k$

${\bf{4}}$ using (4.5) and (4.6) set

and

At each iteration $k$, we may compute the approximate cost

and stop the procedure when

and the desired accuracy (set through an appropriate tolerance parameter $tol$) has been achieved.

Note that the resulting infection dynamics is an approximation of the continuous dynamics governed by the state system (5.6), both because of the optimization algorithm and because of the time discretization of the state equations.

The discretization and the computational procedure outlined above will be used to produce and discuss several scenarios of epidemic control based on the parametrization described in the subsequent section.

7.   Model parametrization and results: application to COVID-19 in Italy

In this section, we parametrize the model using data from the COVID-19 experience in Italy and present the main results. Precisely, we use data from the COVID-19 epidemic from its onset on February 24, 2020 (first data reporting in Northern Italy) until December 19, 2021 when Omicron became dominant, irreparably invalidating the LAR hypothesis. Consistently with our framework, we assume that the course of the epidemic is partitioned into two main periods. The first period (Period 1) covers the initial outbreak (March 2020), the subsequent lockdown (decreed March 20th and ceased in early May 2020), the post-lockdown honeymoon (until early September 2020), and the second pandemic wave, grown during September 2020 and exploding after mid-October due to the collapse of the national tracing system. The second period (Period 2), i.e., the vaccination period, is assumed to begin with the start of the Italian vaccination campaign (January 2021). In previous work ^[17,18], we applied model (5.6) to the COVID-19 epidemic in Italy to identify the different epidemic phases and related parameters. We now apply the algorithm described in Section 6 to compute optimal strategies over the two periods and comparing them with the real course of the epidemic.

7.1.   Data and model parametrization

We hereby summarize the types of data used in the application and the related model parametrization. We report all epidemiological parameters for both Period 1 and Period 2 in Table 1, including vaccine-related parameters (used only in Period 2) and all cost parameters in Table 2.

7.1.1.   Epidemiological parameters

Epidemiological parameters, i.e., those appearing in the state system (5.6), plus those needed to set up the initial data and heredities (5.7)-(5.9), are summarized in Table 1. These parameters will provide the input to (5.6) and (5.11) for the computation of the optimal strategy through the numerical scheme illustrated in Section 6. Most of these parameters are drawn from published work, where model (5.6) was fitted to data ^[17,19].

Note that in Table 1, the value of $I^\#_0$ is drawn from the reported daily cases (see ^[17])

7.1.2.   Cost parameters

The values of cost parameters adopted in the simulations (Table 2) are computed by borrowing input data from a number of published studies and data sources. Given that our model only considers aggregate variables, i.e., it does not include neither chronological age nor other stratifications, a number of intermediate steps were made to compute some of the involved quantities. First, the components of the direct cost $c_E$ of an average infection (see (3.1)) were determined as follows

● the aggregate proportions $g$ of serious infection cases (i.e., cases yielding to hospitalization or death) and the corresponding proportions resulting either in death ($\mu$) or hospitalization without death ($1-\mu$) were drawn from field data on contact tracing collected during the first epoch of the COVID-19 epidemic in Italy ^[27].

● as for the hospitalization cost $\kappa_H$, we considered the aggregate average cost of an hospitalized case (including both the direct cost of hospital stay and subsequent longer-terms treatments after discharge), drawn from a published study on the first pandemic year in Italy ^[28].

● the cost $\kappa_D$ of a death due to COVID-19 was determined by first computing the mean number $T_L$ of life-years lost by an average individual eventually dying at any age with a COVID-19 diagnosis (see Section 3.1). Precisely, $T_L$ was computed by combining age-specific risks of death with a COVID-19 diagnosis from the aforementioned study ^[27] with seroprevalence data collected at the end of the Italian lockdown (^[29]) and with official population data to obtain an age distribution of COVID-19 deaths. The ensuing deaths distribution was combined with age-specific data on life expectancy in Italy to compute $T_L$ as the aggregate (average) number of life-years lost by a person (of any age) dying prematurely due to the epidemic. Finally, $\kappa_D$ was computed by multiplying $T_L$ for the level of per-capita GDP in Italy during 2019, according to (3.1).

As for vaccination costs:

● the average monetary baseline cost of vaccination $\eta$ in equation (3.8) was computed by relying on ^[30], a highly detailed work that estimated the overall per-capita cost of the entire COVID-19 vaccination campaign in New York city by including both the direct costs of the vaccine and its overall indirect costs, including infrastructures, staff, and vaccine administration costs.

● As for parameter $\varepsilon$, since it lacks a direct empirical counterpart, we assume $\varepsilon = 0.001$ as an ansatz.

Finally, as for the indirect costs (function $Q(x)$ in (3.3)):

● parameter $\rho_m$, actually given by the ratio between the level of social contacts in the most altered situation ($c_m$) and its normal level ($c_0$), was approximated by the ratio between the level of transmission under the most intense phase of the Italian lockdown (identified with the lockdown reproduction number ${{\cal R}}_L$) and the transmission during the initial phase of free epidemic (identified with ${{\cal R}}_0$). Specifically, both ${{\cal R}}_L$ and ${{\cal R}}_0$ can be drawn from ^[17] so that, using formula (5.4), we have

● parameter ${\alpha}$ was set to ${\alpha} = 0$;

● parameter $\chi$ was taken as a free simulation parameter;

● parameter ${{\cal L}}$ was assigned based on the loss of GDP observed in Italy during the lockdown period ^[29].

7.2.   Results for Period 1: optimal control in the absence of a vaccine

Here, we investigate the shape of the optimal distancing action over Period 1 only, when no vaccine is available. As a baseline, we continue to stick to the Italian COVID-19 experience and make the following further assumptions on the control horizon:

● the control horizon starts at time $t = 0$ after an epoch of free exponential growth yielding to the number of infective reported in Table 1. The initial date is set as March 15th, 2020, which is a compromise accounting for the progressive completion of the process by which Italy entered in its lockdown, formally declared on March 9th after the closure of school of all grades on March 4th, but further strengthened between March 20 - 25th;

● the length of the horizon ($T_1$) is correspondingly set to $T = T_1 = 307$ days, i.e., the length of the period between the lockdown onset and the start of the Italian vaccination campaign in January 2020.

In particular, we will focus on the dependence of the optimal action on the level of the weight $\chi$ attributed to directs costs. For sake of generality, we will consider $\chi$ values larger or at most equal to the parity level between direct and indirect costs (i.e., $\chi = 0.5$) at which the two cost components are just summed up in the cost functional. Clearly, as the weight attributed to directs costs declines, the intensity of the control action will be attenuated, potentially resulting in a violation of the LAR hypothesis. Cases violating the LAR hypothesis will be discarded a posteriori.

To assist the interpretation of our subsequent results, we first report (Figure 1) the pattern of the overall cost function when a constant control strategy $\rho(t)\equiv \rho$ is adopted, for different $\chi$ values. Notably, overall costs are well-behaved and show a clear minimum point for each $\chi$ considered.

In Figure 2, we compare the optimal strategy identified by the model with the observed trend of the incidence of COVID-19 data resulting from the policy actually implemented by the Italian government during Period 1. In particular, the left panel of Figure 2 reports the temporal courses of the optimal control action $\rho^*(t)$ over Period 1 for different $\chi$ values. These profiles are quite similar with: (ⅰ) a rapid initial decline toward the bottom $\rho_m$ of social activity ("lockdown"); (ⅱ) a maintenance of the lockdown for a duration of about three months (slightly increasing for larger $\chi$ values). The optimal lockdown duration is almost twice as large as the formal duration of the Italian national lockdown (though it should be pinpointed that prudent individuals' behaviors persisted for a much longer time ^[17]); (ⅲ) a rapid increase of social activity to an intermediate level (the same for all $\chi$ values considered), with epidemic reproduction sharply above threshold (with an effective reproduction number ${{\cal R}}_E$ about $1.5$). This intermediate level is maintained for most of the horizon (about five months), unlike the actual Italian policy that abruptly switched from the lockdown to almost normal activity; (ⅳ) a final relapse of social activity (still below its normal level) when the end of the horizon approaches. By comparing the epidemic trends corresponding to the optimal strategies with actual epidemic incidence data consequent to the governmental policy (the dotted curve in Figure 2, right panel), we remark the following main facts: (ⅰ) all the optimal actions resulting from the different $\chi$ values are consistent with the observed data until the half of Period 1 (Summer 2020); (ⅱ) the end of the lockdown phase and the return of transmission above threshold allows infections to restart growing exponentially under all hypotheses considered; (ⅲ) as the end of the horizon approaches and most of the control action is relaxed to keep the indirect costs at a low level, all (predicted) epidemics accelerate sharply. This speed-up of the epidemic at the end of the control horizon is an inevitable phenomenon unless one includes in the cost functional a specific component penalizing epidemic regrowth nearby $T_1$, which we deliberately ruled out; (ⅳ) the most noteworthy effect is that the (seemingly) slightly different dates of exit from the lockdown (resulting from the different values of the weight of direct cost $\chi$) have dramatically different effects afterward. Indeed, at low/intermediate $\chi$ values ($\chi = 0.5, \chi = 0.7$), the return of transmission above threshold brings large epidemics, which become dramatic when the end of the horizon approaches. Instead, at higher levels of prioritization of direct costs ($\chi = 0.90, 0.95$), the –seemingly slight– resulting delays in the optimal exit from the lockdown would allow achieving a substantial containment of the second COVID-19 epidemic wave (debuted in October 2020). In these latter scenarios, the LAR hypothesis is never violated (unless very near the end of the control horizon).

Further insight is obtained (Figure 3) by comparing the optimal social distancing strategy $\rho^*(t)$ for $\chi = 0.90$ with the piecewise, phase-specific, strategy resulting from the multi-phasic approach to the COVID-19 epidemic in Italy proposed in ^[17]. In particular, the (dashed) piece-wise constant curve represents the estimated $\rho$ levels over the different phases of the first years of the COVID-19 epidemic in Italy as identified by ^[17]. Notably, the approach in ^[17] estimated that the harsh control phase initiated with the generalized national lockdown, actually lasted much longer than the official lockdown duration (officially lasted about 50 days and declared over by early May 2020). In ^[17], it was argued that this honeymoon effect was the outcome of a combination of factors including, among other, transmission seasonality, behavioral inertia, i.e., the tendency of people to maintain prudent behaviors even after measures' relaxation, and process inertia. However, with summer time, these protective factors ended and transmission returned much above the threshold, promoting the start of the second epidemic wave (Figure 2, right panel). A comparison of the actually observed level of transmission with the optimal one in Figure 3 is suggestive that the optimal action, implying a more moderate re-opening (30$\%$ less than the actually observed level) and lasting two more months would have allowed fully preventing the second COVID-19 wave in Italy.

7.3.   Including vaccination

We now present the results of the joint optimization of both social distancing and vaccination over the entire horizon $[0, T]$ (given by the union of Period 1 and Period 2) with $T = 644$. The final date of the horizon corresponds to December 15, 2021. By this choice of the horizon duration, we exclude the epoch of the Italian COVID-19 epidemic when Omicron variant became dominant, making the vaccine largely ineffective and therefore causing the LAR hypothesis to be seriously violated. In particular, we will seek the optimal couple $(\rho^*(t), v^*(t))$ minimizing overall costs. Additionally, we will compare the optimal pair with the case in which the observed vaccination schedule is taken as given and set to the actually implemented levels of COVID-19 vaccination in Italy during Period 2. In the latter case, social distancing returns to be the only available control action to minimize costs. In order to make a "fair" comparison between these two cases, we assume that the total number of vaccinations administered over Period 2 in the two cases is the same. This implies setting the bound of the (effective) vaccination rate per unit time to $v_{max} = 0.0029$.

The corresponding results for a high level of prioritization of direct costs ($\chi = 0.95$) are reported in Figure 4. Notably, the optimal vaccination action (right panel) has a full bang-bang shape, i.e., it jumps to its upper bound at the start of the campaign (i.e., at the beginning of Period 2) and terminates at the end of the horizon. This bang-bang shape, which results from the bound imposed on the daily number of administered vaccines, dramatically differs from the complicated vaccination course actually implemented in Italy. The ensuing optimal distancing actions (right panel) shows that, as expected, the fully optimal strategy improves the one based on the observed vaccination schedule, i.e., the corresponding $\rho(t)$ curve lies almost always above the one corresponding to observed vaccination and allows an earlier exit from the lockdown. Additionally, there are several details deserving comments about the optimal social distancing trajectories over the entire period compared to Period 1 (see the left panel of Figure 2). Indeed, before the phase of regular relapse of social activity allowed by the growing vaccination campaign, the optimal strategy over the entire horizon, compared to Period 1 only, emerges as characterized by (ⅰ) a substantially longer lockdown (almost the entire duration of Period 1) and (ⅱ) a dramatically shorter intermediate phase (still about ${{\cal R}}_E = 1.5$) of reduced restrictions. This seemingly inconsistent result is the consequence of redoing a fully independent optimization over the entire horizon by taking the start date of Period 2 as a sure, uncertainty-free quantity determined by the arrival date of the vaccine, which is assumed to be known with certainty. This being said, these differences are amplified by the high degree of prioritization attributed to direct costs. Indeed, by reducing the weight of direct cost to $\chi = 0.7$, both the aforementioned features are sharply mitigated; in particular, the end of the lockdown is backdated by about three weeks (Figure 5), consistently with the findings of Figure 2. Finally, the acceleration in the return to normal social activity nearby the end of the horizon stems from the same arguments reported for Figure 2.

8.   Discussion and concluding remarks

Motivated by the COVID-19 pandemic experience and its momentum on optimal control studies of epidemics, in this work we investigated the optimal control of a threatening epidemic stratified in two main epochs by the availability of control tools, namely a first period (Period 1) during which the epidemic response only relies on NPIs, and a second period (Period 2) during which a combination of NPIs and vaccination is used. With this aim, we proposed a specific numerical algorithm for solving the optimal control of an epidemic model with age of infection developed in our previous work ^[17,18,19]. The cost functional of the model was suitably parametrized by seeking appropriate definitions of the various components of the costs raised by the epidemic (direct, indirect, vaccination) and by combining them with a range of literature data and studies from COVID-19 ^{[27,28,29,30]}. In our analyses, we considered two main problems: a pure problem of optimal social distancing in the absence of a vaccine as well as a problem of joint optimization of social distancing and vaccination under a delayed vaccine arrival. In particular, we focused on the dependence of the optimal solution on the degree of prioritization on direct costs, and compared our optimal solutions with data from the Italian COVID-19 experience.

Our main results are as follows. As for the first problem, the optimal social distancing action shows a stop-and-go behavior, which includes (ⅰ) an initial reduction (bang-bang-like) of social activity to bottom-level (i.e., a "lockdown"), which is maintained for about $30\%$ of the horizon, (ⅱ) then, a partial relaxation of restrictions, bringing social activity to $45\%$ of its normal level lasting for more than half of the horizon, and (ⅲ) finally, a gradual relapse of social activity toward its normal level (but without reaching it), with the purpose of containing the growth of indirect costs as the end of the horizon approaches, thereby allowing epidemic relapse. Comparing our theoretical results with data from the COVID-19 pre-vaccination period, it is suggested that a policy informed by high levels of prioritization of direct costs would have had the potential, by slightly delaying the lockdown termination, to avoid the second (devastating) COVID-19 epidemic wave (debuted in October 2020). As for the full problem, the imposition of a strict bound on the daily vaccination rate forces the optimal vaccination schedule to have a bang-bang shape covering most of the horizon, sharply improving the time-inhomogeneous policy actually implemented in Italy. The corresponding optimal "exit from restrictions" strategy policy is therefore forced to follow a very slow and gradual (about linear) return to normal social activity.

The proposed approach was a parsimonious one that included the age structure of infection in a simple way thanks to a minimal parametrization leading to delay differential equations. Moreover, it relied on the specific assumption that we term the low attack rate (LAR) hypothesis, which showed the ability to describe the COVID-19 epidemic in Italy both before and after the arrival of vaccines ^[17,18] and which can be seen as a way to avoid complications such as the inclusion of saturable public health resources. Further, we proposed a detailed modeling of the different components of epidemics costs and parametrized them by a range of literature estimates. Finally, to solve the ensuing optimal control problem, we proposed an ad-hoc numerical procedure consistent with the theoretical analyses provided in ^[19].

Clearly, the proposed approach suffers a number of limitations, of which we can just list a few here. First, the a priori optimization over the entire control horizon, including a vaccine arriving at a known time, is an obvious drawback that was solved by setting the vaccine arrival time to the observed arrival time of the COVID-19 vaccine. A general solution to this would be to include a random vaccine arrival time as proposed in pioneering economic studies of COVID-19 optimal control ^[7,8]. On the same line of reasoning, the present a priori approach can only be taken as valid for preparedness activities and not for real-time responses to an ongoing outbreak, which would require adaptive control approaches ^[14]. Also, the LAR hypothesis is just a simplification to the true inclusion of the finiteness of public health resources (e.g., hospital beds or testing). The latter issue was first considered long before the COVID-19 pandemic ^[31,32] but has become critical after the COVID-19 onset. Further, COVID-19 has showed a marked seasonality in transmission. The optimal control study of epidemics with transmission seasonality was initiated in ^[33] and should be carefully considered in future preparedness studies. Another issue is that the proposed approach postulates that the decision-maker is able to implement the optimal social distancing action without concerns for individuals' responses. In other words, the present model postulates a fixed population adherence constantly set to $100\%$, which does not need to be the case. Including the responses of individuals to the undertaken policy measures would be an important improvement, which can be carried out by the tools of the modern behavioral epidemiology of infectious diseases ^[34,35,36]. Not to mention about the importance of considering outbreaks evolving due to the onset of variants of the pathogen and the related scenarios of changing vaccine efficacy ^[37].

Notwithstanding the previous limitations, we feel that the present self-contained approach and the ensuing results should be taken as a useful starting point to be improved through the range of ameliorations discussed above.

Acknowledgments

The authors wish to thank the Editors of the special issue "Modelling and computational strategies for infectious disease dynamics and ecology of structured populations" for considering the present work in the volume, and the referees for their suggestions that improved the manuscript.

MI and GM gratefully acknowledge the support received in the framework of the Italian-Romanian collaboration agreement "Analysis and control of mathematical models for the evolution of epidemics, tumors and phase field processes" between the Italian CNR and the Romanian Academy.

G.M. acknowledges the support of a grant of the Ministry of Research, Innovation and Digitization, CNCS - UEFISCDI, project number PN-Ⅲ-P4-PCE-2021- 0921, within PNCD ${\rm{IIII}}$.

Use of AI tools declaration

The authors declare that they have not used Artificial Intelligence tools in the creation of this article.

Conflict of interest

All authors declare that there are no competing interests.