Long-term spatial and population-structured planning of non-pharmaceutical interventions to epidemic outbreaks.

In this paper, we consider the problem of planning non-pharmaceutical interventions to control the spread of infectious diseases. We propose a new model derived from classical compartmental models; however, we model spatial and population-structure heterogeneity of population mixing. The resulting model is a large-scale non-linear and non-convex optimisation problem. In order to solve it, we apply a special variant of covariance matrix adaptation evolution strategy. We show that results obtained for three different objectives are better than natural heuristics and, moreover, that the introduction of an individual's mobility to the model is significant for the quality of the decisions. We apply our approach to a six-compartmental model with detailed Poland and COVID-19 disease data. The obtained results are non-trivialand sometimes unexpected; therefore, we believe that our model could be applied to support policy-makers in fighting diseases at the long-term decision-making level.

Introduction

In the face of epidemic respiratory viruses, when medical resources are not yet available, the very first response methods to delay and moderate their spread in the population are usually concentrated on non-pharmaceutical interventions (NPIs) (Haug et al., 2020). NPIs include interventions, such as contact tracking and social distancing, including stay-at-home recommendations, the closure of leisure activities, schools and day-care centres, and the wearing of face masks, as was the case with COVID-19 in 2020 (Djidjou-Demasse et al., 2020, Yang et al., 2021). Given the high economic and social costs of NPIs, it is critical to understand whether these interventions have the desired effect of controlling the epidemic, and which interventions are necessary to maintain control. Before vaccines were introduced, NPIs were the only measure to slow down the spread of infection. After the introduction of vaccines, and despite their availability, maintaining NPIs was still critical in limiting the consistently high number of infections (Moore et al., 2021). The effect of major interventions on the time-varying reproduction number (Rt) was studied across 11 European countries for the period between February 2020 and 4 May 2020. The results suggest that, taken together, these interventions have had a substantial effect on transmission, as measured by changes in the estimated Rt. It was found that estimates of Rt range from a posterior mean of 0.44 for Norway to a posterior mean of 0.82 for Belgium, with an average of 0.66 across the 11 countries, which account for an 82% reduction compared to pre-intervention values. For all countries, it was found that NPIs have been sufficient to drive Rt below 1 and achieve control of the epidemic (Flaxman et al., 2020). The obtained results correlate with studies from individual countries like Spain (Hyafil and na, 2021), France (Salje et al., 2020), Germany (Robert Koch Institute, 2020) and the UK (Davies et al., 2020) conducted over a similar period, which arrive at very similar estimates despite different methodologies and data. The NPIs had a very significant estimated effect on deaths. By comparing the deaths predicted under the Seth Flaxman et al. model with no interventions to the deaths predicted as a result of the intervention model, the total deaths averted over the study period was calculated (Flaxman et al., 2020). It was found that across 11 countries, about 3.1 million deaths have been averted owing to the introduction of interventions.

Previous outbreaks along with relevant research have shown that compartmental epidemiological models help to understand and make decisions on the response to the spread of disease (Holmdahl and Buckee, 2020). The SIR model captures the basic epidemiological features of infectious diseases, that is the transition of the population from susceptible (S) to infectious (I) and finally recovered (R) states along with the relevant contact- and disease-specific rates (Kermack et al., 1927). In order to enhance quality, more complex, extended SIR-type models are developed and analysed, challenging in particular the assumption on homogeneous mixing  (Xiang et al., 2021).

Even though epidemic models proved to be useful in the short term to plan interventions or to investigate their impact, yet their applications are limited in a longer time horizon because they do not incorporate feedback on how the pandemic spread may change. That is why a stream of research, which includes this paper, is devoted to implementing control decisions directly into the epidemic model. Instead of focusing on predicting the future epidemic process, optimisation models determine disease response strategies in line with the direction in which epidemic processes should evolve. This is especially important for decision makers, who are constrained and want to apply the least restrictive programmeme (Fenichel, 2013), or who have limited access to resources such as vaccines etc. (Dasaklis et al., 2012).

Constrained decisions can only serve as the desired trade-off solution between economics and infection dynamics when they are targeted (Glass et al., 2006), which means that modelling should be detailed, keeping the computational complexity in mind. Spatial and demographic heterogeneity models are believed to bear the closest relation between targeted interventions and the dynamics of any epidemic (Lloyd and May, 1996, Greenhalgh, 2013), and particularly COVID-19 (Badr et al., 2020, Chinazzi et al., 2020, Khajji et al., 2021, Medrek and Pastuszak, 2021). The aim of our work is to develop a model that considers the mobility and structure of regional populations when optimising non-pharmaceutical interventions, that restrict transmissivity between susceptible and infectious individuals. We strove to show that such a model enables us to make better targeted decisions, and can be used in reasonably large scale cases. In spite of their potential usefulness and importance, practical models of this kind can be scarcely found in the literature, because the detailed representations of disease transmission processes can make the optimisation of controls too complex (Bussell et al., 2019).

Contributions. We formulate a new model that optimises strategies to control the spread of epidemic disease by efficiently planning non-pharmaceutical interventions in the long term. The main contributions are:

We introduced optimisation decisions into a discretised version of the compartmental epidemic model that considers the detailed network of contacts between age-structured groups located in different regions. According to our model, infection in a specific location develops due to contacts between residents and visitors from other regions. Specifically, we consider the possibility of new infections occurring due to contacting individuals from other areas. Thanks to this, we can deliver a targeted, multiperiod plan for interventions that accurately impact mobility to contain the epidemic.

The results we provide prove that the detailed modelling of contacts in optimisation models is justified and necessary, and that the quality of planning is increasing along with spatial model granularity. This brings a computational burden. We show that the resulting large-scale non-linear and non-convex model can be solved with a specific variant of covariance matrix adaptation evolution strategy and gives better solutions than heuristics.

A combination of the two previous aspects, model formulation and solution methodology, allows us to investigate the impact of spatial model granularity, on the one hand, and the role of policy makers’ preferences reflected in the three variants of the objective function on the other. The real data for Poland in case of COVID-19 were used to justify our approach, showing that the Polish model with 380 counties offers the most reasonable trade-off between the quality of planning and computational complexity. We show that computed policies can lead to different results, depending on criteria; for instance, minimising the maximal peak in new infections may not result in a minimal number of deaths. No other works, to the best of our knowledge, provide such a detailed insight into the planning of action strategies for epidemics on such a large scale nationwide. That is the reason we believe that our model could play an important role in a wider spectrum of OR tools supporting policy-makers during epidemic outbreaks.

The paper is organised as follows. Section 2 reviews the related literature with special consideration to spatial epidemic models and disease control models, which are rare. In Section 3, we provide a detailed description of a spatially heterogeneous SIR-based model used as a basis for developing constraints in an optimisation model, as explained in Section 4. More specifically, we explain how to transform the formulas to define control variables for each area and age group, while capturing a spatially heterogeneous mix of individuals at the same time. We also discuss the objective functions and constraints related to limits on restrictions and resources, and the solution algorithm. In Section 5, we apply our method to analyse spatially heterogeneous strategies for the COVID-19 outbreak in Poland. Finally, we discuss the robustness of the solution in Section 6. Concluding remarks are given in Section 7.

Literature review

The need to use network classes to model populations in order to study epidemic processes is thoroughly explained in a survey by Pastor-Satorras et al. (2015). They review classical epidemic models and the characteristics of complex networks, and discuss how different methodologies address heterogeneous contact patterns. The authors concentrate on approaches considering high granularity networks, in which each node corresponds to a single individual linked to others in a population. One stream of research originates in cyber-security problems, and attempts to recognise links in order to remove them from the network with infected nodes. In an epidemic context, Nandi and Medal (2015) propose mixed integer linear programme models representing connectivity-based networks with infected and susceptible nodes. If the highest probability transmission paths are removed, then the infection speed is reduced, which is shown on randomly generated networks, each with 12 nodes. Individual contact networks are especially useful to investigate immunisation strategies, see for example the work by Glass et al. (2006). At the level of the individual, the age of the person may be easily taken into account, and that is one of the most important characteristics related to the effectiveness of decisions that can mitigate epidemics, as surveyed by Greenhalgh (2013) in age-structured models.

The first attempts to incorporate metapopulation mobility processes into the SIR epidemic model were proposed in the 90s (Grenfell and Harwood, 1997, Lloyd and May, 1996, Sattenspiel and Dietz, 1995). Researchers noticed that contact between individuals occurs as a result of the mobility of participants across either geographic or social space. They divided a population into subgroups homogeneous in their mobility patterns (geographic, social mobility etc.). The force of infection depends on the transmission from locals (both travelling and non-travelling), as well as from visitors. Our approach is similar, which is rare in models optimising epidemic spread.

Recently, several researchers have developed analytic and forecasting models to assess the impact of a number of control strategy scenarios.  Mistry et al. (2021) developed age-stratified contact matrices for 35 countries, and used them to model the spread of airborne infectious diseases. Their results show that sub-national heterogeneities in human mixing patterns have a marked impact on epidemic indicators. Interestingly, as Chinazzi et al. (2020) have shown, when mobility is aggregated modelling results indicate that sustained 90% travel restrictions to and from mainland China only modestly affect the epidemic trajectory unless combined with a 50% or higher reduction of transmission in the community. The results of Gatto et al. (2020) suggest that the sequence of restrictions posed to mobility and human-to-human interactions have reduced transmission by 45% (42% to 49%). Balcan et al. (2009) used a gravity law to model the number of visitors travelling across geographical locations and developed a hierarchical framework considering both short-range mobility and long-distance travels, showing that both levels should be carefully addressed by different mobility restrictions. At the global level, the modelling of a disease spread over time and spatial scales is proposed by Balcan et al. (2010), where the age-specific contact matrix reported by Wallinga et al. (2006) was incorporated into computational simulations. Other approaches may be found, see for example Badr et al., 2020, Giordano et al., 2020, Medrek and Pastuszak, 2021 and Pei et al. (2020). A growing amount of research confirms and highlights the need to address mobility and age groups when evaluating interventions. It is difficult, however, to include them in optimisation models due to their high dimensionality and the large scale of epidemic models.

A simplified approach to mobility is to assume that the population size in each region remains constant for every time step, while transmission is relative to the proportion of contacts between each pair of regions. Such an approach was proposed by Zakary et al., 2017a, Zakary et al., 2017b to derive optimum control strategies based on a direct definition of the necessary conditions, which leads to the calculation of the Hamiltonian derivative. These works carried on in the COVID-19 context, i.e. by Hamza et al. (2020) or Moyo et al. (2020). The underlying methodology is hard to be applied on a large scale. Only small cases (3 to 5 regions) were solved.

Studies that combine spatial and population-structured optimisation models with dynamics of epidemic outbreaks more often appear in the operation research stream devoted to medical resource allocation. These typically include vaccine allocations, but under the recent COVID-19 pandemic a wide variety of logistics operations are examined, i.e. personal protective equipment allocation (Dönmez et al., 2022), ventilator distribution (Yin et al., 2021), or patient allocation to facilities (Hosseini-Motlagh et al., 2021); see also the review by Queiroz et al. (2020). These studies derive discrete-time epidemic models to include them in resource allocation models. This paper falls into the same limited area of studies that implicitly consider a discrete epidemiological SIR-type metapopulation spatial model within the optimisation problem.

One recent approach motivated by optimising vaccine and antibiotic drug allocation decisions for a cholera outbreak was proposed by Du et al. (2021). The authors verified their three-step data-driven approach on Haitian data divided into ten departments. They compared the optimisation-based allocation strategy’s humanitarian and economic performance with two commonly used strategies: static and reactive. Another recent research in this field is the work by Abdin et al. (2021) who developed a novel epidemiological model for allocating limited testing capacities within different segments of the population. A non-linear programming model is solved using the Interior Point Optimiser (IPOPT) solver for three interconnected regions in France under the COVID-19 pandemic outbreak. A wide variety of scenarios was analysed to provide insights on optimal timing, location and level of tests to be distributed. Finally, Rezapour et al. (2021) developed a bi-level optimisation approach, where the strategic decisions of determining the least costly mobility restrictions for a metropolis are imposed onto a so-called reaction–diffusion process developed to model infection transmission. A model minimising socioeconomic costs is tested and compared with two intuitive intervention schemes in the Sioux Falls metropolitan area in the U.S. (24 nodes and 38 links).

Our approach answers a novel research question within the prior literature. What benefits can be brought from a more detailed network model of disease spread when optimising non-pharmaceutical interventions. We strive to indicate a good trade-off between optimisation capabilities and the quality of obtained decisions. To the best of our knowledge, none of the previous studies investigated the combined NPI planning model and SIR-based metapopulation spatial model on a large-scale, nation-wide level.

Spatial SIR-based model

Spatial SIR-based epidemiological model

Building a decision model for NPIs requires the incorporation of an epidemiological model that translates decisions into the observed effects of these decisions. Markov chain-like compartmental models are widely used to model the progress of infectious diseases. The SIR model is the best-known, most basic model in epidemiology (Kermack et al., 1927). It assumes three compartments: —individuals not yet infected, but susceptible to the disease at time , —individuals who have been infected and are capable of spreading the disease at time , and —individuals who have been infected, but at time are removed either due to immunisation or death. The dynamics of transition among compartments is given with differential equations: where is the transmission rate per contact, and is the recovery rate. There are several SIR-derived models in the literature which differ in number and the definitions of compartments; however, they usually assume full homogeneous mixing of individuals in the population and constant population , as the SIR model does. Neglecting the demographic transition of birth and death is justified for a fast-spreading disease, the homogeneous full mixing assumption may significantly disturb the results. Therefore, we build our disease spreading model based on the SIR-based family, but we incorporate into the model the heterogeneous mixing of individuals. We assume that the mixing of individuals depends on their spatial interconnections, as well as existing subpopulations; however, in specific locations mixing is homogeneous (see Fig. 1).

Fig. 1

Model of contacts considered in this paper takes into account spatial heterogeneity of subpopulations, e.g. age groups, while mix between infectious and susceptible people in a specific location is homogeneous.

Spatial interactions are mostly caused by daily mobility associated with travel to workplaces and schools, which usually constitutes the main steady stream of mobility. Taking that fact into account, let us model the spatial interactions of individuals by the directed graph , where represents the spatial division of space into regions where individuals exist and represents possible interactions between individuals from two regions connected by edge . We assume discrete time, typically days; therefore, the considered horizon is divided into time intervals. The mobility factor is associated with each edge . It is a fraction of the population of the source vertex, travelling to the destination vertex along vertex at the beginning of each time interval. We assume that these individuals spend their entire time interval at their destination location and at the end of it they travel back to their source. We neglect travel time and interactions that may occur during travel. Obviously, for some edges, the mobility factor can be equal to zero; therefore, the resulting graph does not have to be a full graph.

We assume full mixing of individuals at each node within a time interval. In particular, the mixing of individuals available at the beginning of interval results in new infections occurring at the end of the interval, i.e. . We do not neglect the fact, however, that during each time interval at vertex there are individuals coming from different locations, mixing all together. This includes the possibility of infectious individuals travelling to a destination and spreading the disease to individuals not only originated from that node, but also onto individuals from different locations visiting this node during this time interval. Some papers propose the application of an SIR-based model for each node separately, but they simplify the interaction model by assuming a constant population at each node, as the original SIR model does. We model these details and take into account the fact that the temporal population at node varies, since it depends on the number of visitors. The number of individuals who transforms from compartment to is the sum of transformed individuals from over visiting destinations. Therefore, under the assumption of discrete time , where is a set of periods in the horizon under consideration, we transform model (1)–(3) into a differential one: where is the fraction of individuals from visiting over interval , is the effective number of individuals at and over interval , is the effective number of infected individuals at and over interval , while are the numbers of individuals in each compartment at at the beginning of interval . For pairs of nodes between which there are no travellers, the coefficient takes the value of 0, so that the above equations can be stated in a general form. In practice, matrix can be quite sparse.

The effective number of infected individuals is the result of visiting infected individuals as well as local infected individuals, taking into account that some of them visit other nodes during interval . The number of infected individuals can be calculated as follows: Since the sum also includes the factor, is equal to the total number of infectious individuals. Similarly, can be calculated in the following way: where is the constant number of individuals at at the beginning of each time period . For all , , since individuals are travelling to other destinations for period , but they go back at the end of each interval .

We also consider population branching, that is, we assume that a population is divided into subpopulations, e.g. age groups. Let be the set of subpopulations. The mobility factor is defined for each subpopulation and denoted by . This allows us to reflect differences in the mobility of each subpopulation. For instance, some important mobility patterns can be expected when age groups are considered, since working people travel differently than students or seniors. Moreover, we allow for edges , which means that mobility factor is a fraction of subpopulation that mixes with other individuals at their origin vertex . This allows us to model the fact that not every individual mixes with others on a daily basis, especially taking into account population branching.

Population branching requires the reformulation of (4)–(6) to subpopulation-specific equations. From the perspective of individuals mixing, is the key element in the model. Let us define as a number of individuals transforming from compartment to compartment at location , subpopulation and over time interval , defined as follows: Then, the spreading model can be stated as follows:

Six-compartment spatial SIR-based model

The formulation we proposed in the previous subsection can be applied for any SIR-derived model. It is known that there are infectious diseases with relatively long incubation time. This means that there is a significant latency period after infection with disease and becoming infectious. After that period an individual may become symptomatic or may go through the disease asymptomatically. Therefore, let us consider a class of diseases with significant incubation period and symptomatic/asymptomatic cases. COVID-19 is an example of a disease that suits this class well. To reflect the properties of such a disease let us apply our approach to six-compartment model. To take into account the relatively long latency period, we introduce an exposed individuals compartment (E). We also introduce an asymptomatic patients compartment (A), while compartment I remains for symptomatic cases. Additionally, the model differentiates compartment R into recovered (R) and dead (D). This type of model, denoted as the SEAIRD model, has been proven to accommodate accurately the actual data of the COVID-19 epidemic in various countries (Giordano et al., 2020, Liu et al., 2021, Roy et al., 2021, Xiang et al., 2021). Finally, the spatially heterogeneous transmission SEAIRD model can be stated as follows: where is incubation rate (inversion of incubation time), is fraction of exposed individuals who become symptomatic, and are recovery rate in case of symptomatic/asymptomatic individuals respectively, is mortality rate at subpopulation , and are respective transmission rates per contact of symptomatic and asymptomatic individuals. The fraction of symptomatically infected individuals from subpopulation living at node , who travel to node is denoted by . The parameter was introduced to distinguish the travel rates of symptomatic and asymptomatic individuals. The population size at node that takes part in the mixing is computed taking into account all the traffic between the nodes. It is computed as the total number of all individuals travelling to including those already living in .

Optimisation model

The model introduced in the previous section lays the background for the decision model, which is the primary goal of this paper. Decisions are related to non-pharmaceutical interventions that may include for instance the mandatory wearing of masks, limits in the number of persons in public spaces, social distancing requirements and the closure of selected facilities. In general, we assume that such decisions can be made for each node and subpopulation separately, although in the next section we consider some aggregation of decisions. Let us introduce which reflects the level of non-pharmaceutical interventions at location applied to subpopulation in an interval . Intervals used for decisions can aggregate intervals from set used in the disease spread model. Let be mapping , which maps time interval into interval . Typically, while the spread of disease is expected to be modelled for daily granularity, decisions can be modelled weekly, since it would be hard to apply any policy that changes the restrictions more frequently.

In general , where means no interventions and full separation of individuals. In reality, must be bounded at every period as follows

Since reflects a kind of oppressiveness, we limit the average value of as follows:

Any reduces transmissivity of disease during contact. Therefore, we include our decision variables in (19) in the following way:

The resulting optimisation model takes the following form: subject to where objective function reflects the required optimisation criteria. We consider three variants of objective function:

sumI—total number of symptomatic cases, which reflects the natural willingness to limit the fraction of the affected population:

maxI—maximum number of symptomatic cases, which is related to the required resources in hospitals, like the number of beds or respirators. This requires the introduction of an additional variable , which is the maximum of cases over time which should be minimised:

sumD—total number of deaths directly resulting from the development of the disease:

Each variant of the model is a non-linear non-convex continuous optimisation problem. Moreover, as we show in the next section, it is important to include high spatial granularity of the spreading model and decisions, so that the model becomes a large-scale optimisation task. The case of Poland that we discuss in the following section assumes three group ages, a six-month horizon and 380 regions, which gives more than 207 thousand primary variables . According to the data we possess, the spread model, as well as decisions, could be potentially modelled for a more detailed administrative grid, but this would result in a tremendous number of variables, namely, 1352 thousand primary variables.

High non-linearity and the large number of variables do not allow for the direct solution of the mathematical programme. We treat the spread model (13)–(18), (23) as a blackbox function and we introduce constraint (22) to the objective as a quadratic penalty function, which results in the following unconstrained large-scale optimisation problem: where is the penalty weight. Although constraint (22) is considered as a soft constraint, the numerical experiments presented in the following sections show that it was satisfied near equality.

In order to solve such a reformulated problem we applied a variant of covariance matrix adaptation evolution strategy (CMA-ES), which is considered as state-of-the-art in evolutionary computation for large-scale problems (Hansen and Ostermeier, 2001). The algorithm is known as a separable CMA-ES, and it limits the covariance error update to the diagonal elements, which reduces precision at the benefits from lower complexity that becomes linear in the variables space dimension. In terms of the original algorithm, it is quadratic (Ros and Hansen, 2008). More specifically, we applied the libcmaes library developed at the Laboratory for Computer Science at Université Paris-Sud (https://github.com/CMA-ES/libcmaes), authored by Emmanuel Benazera and supported by the coauthor of the original method, Nikolaus Hansen. All solutions were obtained with the size of the solution population set automatically by the library. The starting solution was set to no interventions, that is, . We used three termination criteria: (1) objective tolerance (1E−7), (2) maximum number of objective evaluations (2E6 unless otherwise stated), and (3) computation time limit (8 h). Any time the algorithm finished with either criterion (2) or criterion (3), we emphasis it in bold, throughout the paper.

Computational experiments

Spread model variants

We consider four different spatial models based on the administrative divisions of Poland:

Country—one node representing the whole of Poland;

Voivodeship—16 node model, each node representing a province (known as voivodeships);

County—380 node model representing counties;

Community—2477 node model representing communities (municipalities);

The number of inhabitants in each age group are taken from Poland Statistics reports. Mobility between nodes is estimated based on mobility statistics for employees and students that are known for the Community model. Based on that, they are aggregated for other models (Statistics Poland, 2018, Statistics Poland, 2019). Throughout the paper we assume that the Community model is the one that reflects reality in the best way. Any other model is a result of aggregation of data available in the Community model.

The initial state is defined by the number of cases in each compartment at time . We consider two initial state variants. The first one assumes that distribution is even across all nodes. Quite often, we deal only with the aggregated number of cases for the entire region and assume even distribution. In the second variant, we calculated the number of symptomatic cases in each community in relation to the reported cases from one particular day. Then, we calculated the number of individuals in each compartment and at each node according to the global proportion to symptomatic cases. This gives us an uneven distribution, reflecting the number of reported cases from the referent day.

The complete list of values of parameters is presented in Table 1.

Table 1

Parameter values.

Six-compartment SIR model (SEIARD) parameters
Param.	Definition	Value	Reference
β	Per contact rate of transmission from infected symptomatic	0.2	Medrek and Pastuszak (2021)
μ	Per contact rate of transmission from infected asymptomatic	0.2	Medrek and Pastuszak (2021)
γ_I	Recovery rate of infected with symptoms	1/21 (0.04762)	WHO (2020)
γ_A	Recovery rate of asymptomatic individuals	1/14 (0.07143)	WHO (2020)
σ	Incubation rate	1/5 (0.20)	WHO (2020)
ρ	Fraction of individuals who become symptomatic	0.2	Giordano et al. (2020)
η^Y	Mortality ratio in young group	7.2E−7%	WHO (2020) and Statista (2020)
η^M	Mortality ratio in middle-aged group	1.5E−4%	WHO (2020) and Statista (2020)
η^O	Mortality ratio in elderly group	7.3E−3%	WHO (2020) and Statista (2020)
Initial condition parameters given in relation to the size of reference population N
Param.	Definition	Value	Source
I(0)	Percentage of infected population	6.5E3/(γ_I⋅N)	Polish government statistics
A(0)	Percentage of asymptomatic population	40E3/(γ_A⋅N)
E(0)	Percentage of exposed population	50E3⋅σ/N	Available at https://www.gov.pl/
D(0)	Percentage of deaths	43.7E3⋅/N
R(0)	Percentage of recovered	4.2E6/N

Heuristic strategies

Constraint (22) limits the total NPIs and, therefore, it can be said that it defines a certain type of NPI “budget”, which may be used in heuristic algorithms. As a reference point for optimised strategies, we define five heuristic strategies as follows:

noNPI—no non-pharmaceutical interventions applied

Evenly—the whole “budget” of NPIs is distributed evenly over time, age groups and spatially

Greedy—NPIs are set to maximum until the “budget” is exhausted

DeathRatio—NPIs are distributed among age groups proportionally to the death ratio, evenly across space, and greedy with respect to time

GreedySpatial—NPIs are set to maximum for regions with the highest infected ratio until the “budget” is exhausted

Heuristic noNPI is only used as a reference point to any other strategies that apply NPI and is expected to be always the worst one. The Evenly and Greedy heuristics are simple strategies that do not differentiate decisions spatially. DeathRatio tries to check whether differentiation among age groups brings meaningful benefits. GreedySpatial is a strategy that uses knowledge about cases to allocate the largest interventions to regions that are highly infected. Aggregated NPI controls are illustrated in Fig. 2. Since the DeathRatio heuristic is the only one that differentiates controls in group ages, it is presented per group in the second chart.

Fig. 2

Aggregated controls for heuristics in the Community model. Second figure shows controls for DeathRatio strategy for each age group.

For each strategy, we assumed a horizon of 182 days, which is reasonable for policy-makers. It is long enough to be a valuable source for making plans according to business and industry. At the same time, any longer horizon would not be credible, bearing in mind the great variety of factors that may shape the development of the disease.

The number of symptomatic cases for different strategies and model granularity are presented further in Fig. 4. Fig. 3 presents all model compartments for the Community model.

Fig. 4

Number of symptomatic cases for different granularity and heuristics.

Fig. 3

Compartments I, A, D, E and S, R of the Community model under no NPIs.

Impact of granularity of disease spread model

Bearing in mind that the Community model is our best approximation of reality, we simulated the disease spread with models of lower granularity and compared the obtained results to the Community model. Table 2 presents deviations from the simulation results with Community for three measures, e.g. sumI, maxI, sumD, under the assumption that initial disease distribution, understood as a number of cases at time 0, is even. For any spatially uniform strategy, the result does not depend on granularity, and in such a case there is no advantage in using a more detailed model. This led us to the following proposition.

Table 2

Total and maximal number of infected, and the total number of deaths for even initial distribution. Percentages show deviation from the Community model.

	noNPI
	Country	Voivodeship	County	Community
sumI	0.00%	0.00%	0.00%	6 239 752.91
maxI	0.00%	0.00%	0.00%	1 816 813.82
sumD	0.00%	0.00%	0.00%	128 796.61
	Evenly
	Country	Voivodeship	County	Community
sumI	0.00%	0.00%	0.00%	5 348 720.21
maxI	0.00%	0.00%	0.00%	1 113 366,61
sumD	0.00%	0.00%	0.00%	109 140.76
	Greedy
	Country	Voivodeship	County	Community
sumI	0.00%	0.00%	0.00%	5 630 358.89
maxI	0.00%	0.00%	0.00%	1 143 746.53
sumD	0.00%	0.00%	0.00%	105 926.97
	DeathRatio
	Country	Voivodeship	County	Community
sumI	0.02%	0.01%	0.01%	5 423 970.83
maxI	0.09%	0.06%	0.02%	927 063.10
sumD	0.16%	0.12%	0.05%	95 917.76
	GreedySpatial
	Country	Voivodeship	County	Community
sumI	13.92%	−5.47%	−1.08%	4 942 460.48
maxI	0.53%	−4.04%	−1.16%	1 137 661.62
sumD	4.17%	−6.66%	−0.72%	101 694.66

Let us consider model (13) – (18) , (23) with constant initial ratios of each component. Any strategy that is spatially uniform leads to the same spreading process for any granularity of the spread model independently of the intensity of an individual’s mobility.

We provide a sketch of the proof. Let us consider two disconnected nodes. Since the mix of compartments and the spread ratio are the same at each node, the disease spreads at the same pace at each node. If we assume full mixing of all individuals from two nodes, then the mix of compartments and the spread parameters stay the same, so the pace of spread is not affected and the results do not change. Now, let us assume that a number of individuals travels from the first node to the next. Since the ratio of asymptomatic individuals in the travelling group is the same as at their origin and destination, this does not change the pace of spreading in any nodes, and the development of the disease is not sensitive to the shifts of individuals. □

While strategy DeathRatio reveals a relatively small impact of differentiation on age groups, strategy GreedySpatial shows that when we allow for a differentiation strategy across nodes in the model, the impact is significant, resulting in nearly 14% error in relation to the Community model.

Table 3 presents similar results, but now the initial distribution of compartments is uneven. In that case, differences also emerge in spatially uniform strategies and, in general, are bigger. An interesting fact is that deviations can be negative or positive. In order to explain that, first let us introduce the following proposition.

Table 3

Total and maximal number of infected, and total number of deaths for uneven initial distribution. Percentages show deviation from the Community model.

	noNPI
	Country	Voivodeship	County	Community
sumI	−0.33%	−0.34%	−0.34%	6 260 269.60
maxI	12.42%	11.02%	8.28%	1 616 081.32
sumD	0.41%	0.43%	0.37%	128 275.50
	Evenly
	Country	Voivodeship	County	Community
sumI	0.21%	0.17%	0.19%	5 337 777.94
maxI	13.10%	11.20%	7.99%	984 390.24
sumD	3.73%	3.68%	3.35%	105 226.80
	Greedy
	Country	Voivodeship	County	Community
sumI	−1.34%	−1.45%	−1.27%	5 706 614.94
maxI	−21.76%	−19.77%	−16.06%	1 461 773.92
sumD	4.20%	4.17%	3.86%	101 666.72
	DeathRatio
	Country	Voivodeship	County	Community
sumI	−2.69%	−2.59%	−2.17%	5 575 062.88
maxI	−26.80%	−24.45%	−20.03%	1 267 678.58
sumD	3.44%	3.52%	3.30%	92 874.16
	GreedySpatial
	Country	Voivodeship	County	Community
sumI	8.52%	−7.53%	−2.96%	5 188 250.84
maxI	−5.68%	−9.13%	0.07%	1 212 636.53
sumD	1.74%	−4.68%	−0.52%	104 128.57

Let us consider the network SIR model with two nodes, A and B. The number of susceptible and infected individuals are denoted by and, respectively. For simplicity, we assume that the number of removed is 0 at each node. We consider two cases: disconnected separated nodes and full mixing. The pace of disease spread in the case of full mixing is not lower than in the case of separated nodes.

In the case of disconnected nodes, the pace of disease spread is . In the case of full mixing, the pace is . Then, which is always non-negative, so . □

Fig. 4 shows the number of symptomatic cases for different granularity. According to Proposition 2, we observe more intense disease spreading, which results in a more sloped function . More intense course of the disease, however, also means faster disease suppression. Assuming a limited horizon of simulation, the relation between these two factors results in deviation that can be either positive or negative.

More granular models better approximate the results achieved with the Community model; however, the gap stays significant. NPI policy optimisation should be applied for the most detailed network model, especially if the initial spread of disease is uneven and non-uniform policies are taken into consideration. Conversely, there is an efficiency barrier that does not allow us to run optimisation on very detailed models. The question is whether it is better to use a less detailed network model and optimise decisions better or to use a more detailed network model at the cost of a worse optimisation process. Therefore, there is a need to search for a good trade-off between the spatial granularity of disease spread models and the capabilities of optimisation techniques. What is also extremely important, however, is that our results show that any further comparison of NPI policies needs to be conducted on the same granularity models.

Impact of decision granularity

In the previous section we showed that network granularity is meaningful. Optimisation on a very detailed spreading model could be possible if the decisions were less granular; however, the question will be the impact on optimised NPI policy. We ran optimisation for Country, Voivodeship and County models assuming that decisions are optimised with the same granularity as the disease spread model, that is for the whole country, for voivodeships and for counties, respectively. Then, we compared it to results of optimisation under the assumption that decisions are always aggregated for the entire country. The results are presented in Table 4. Of course both cases should yield the same results in the case of the Country spreading model, which means that spreading and decision granularity reflect a one-node model. For the Voivodeship model, difference between disease spreading granularity at voivodeship level and spatially uniform decisions is around 5%–8%, while in the case of aggregation decisions from the county level to spatially uniform level is around 2%–4%. Therefore, we can conclude that the aggregation of decisions may worsen the results by several percent. What is interesting is that the biggest impact has been observed for the sumD criterion.

Table 4

Total number of symptomatic cases, deaths and maximum of daily symptomatic cases for spatial optimisation and relative deterioration when spatially uniform decisions are optimised.

	Country		Voivodeship		County
sumI	4 757 137.31	0.00%	4 570 455.31	5.04%	4 729 246.61	2.73%
maxI	594 674.52	0.00%	557 364.86	6.09%	580 260.63	1.80%
sumD	88 993.73	0.00%	82 468.87	8.23%	85 732.94	4.33%

Fig. 5 shows the distribution of new cases and deaths aggregated over the entire time horizon. For the objective sumI, most cases are found in the group of middle-aged individuals, which is justified by the relatively higher mobility between voivodeships (trips to schools are rather local). In the case of sumD, deaths are mainly spatially spread over the elderly group of individuals. At first glance, the solutions seem to be highly correlated, almost proportional. It can be observed, however, that relative differences of objectives are not very different, but absolute differences are still significant. The total number of deaths is equal to 86,520 in the case of sumI, compared to 82,468 in the case of sumD, while the total number of symptomatic cases is 4,570,445 versus 4,621,276 in the sumI and sumD cases respectively.

Fig. 5

Distribution of new cases and deaths at Voivodeships for sumI and sumD criteria.

Improvement possibilities in relation to heuristics

In order to understand what one may expect from optimised NPI strategies, we compared our optimisation results with simple heuristics. Based on the conclusions from previous sections we assumed that decisions are made on the same granularity as the disease spread model and we compared the results with the Community model as our best approximation of the real world. In particular, the Country model means that the whole country is modelled as one node and decisions are uniform all across the country. Policy obtained from the optimisation process is used for a simulation with Community granularity. For the County, model we applied the disease spread and decision models at county level, and the result of the optimisation process was used in asimulation at the community level.

NPI heuristic policies are also simulated at the community level. For instance, the Greedy strategy for the County model means the allocation of NPIs according to the Greedy algorithm applied to the county level, and then the application of that strategy to the Community simulation model. For each optimisation criterion, we choose the best result that we obtained by all heuristics and that value is used as a reference for optimised results. This approach guarantees that obtained results are comparable and verified against our best approximation of the real spreading process.

Table 5 presents the results of the comparison. In every case we managed to improve the results. The best results are for the County case, which may mean that the potential for optimisation in such a model is the largest. For instance, 18.17% in terms of sumD means saving 16.9 thousand individuals, while 23.11% in terms of maxI means that disease peak is lower by almost 249 thousand cases, which highly reduces the required number of beds in hospitals.

Table 5

Comparison of optimisation results against best results obtained by heuristics.

	Country		Voivodeship		County
sumI	5 222 375.68	6.33%	4 434 413.12	0.60%	4 727 407.56	1.99%
maxI	1 095 608.92	4.21%	829 756.09	8.58%	828 225.91	23.11%
sumD	87 531.64	5.75%	81 540.07	10.37%	75 997.53	18.17%

Fig. 6 illustrates how the aggregated NPI control strategies change for different objectives. It can be noticed that for the maxI objective the strategy is to introduce high restrictions quickly and then gradually reduce them by the end of the considered time horizon. This results in a number of increasing symptomatic cases at the beginning, but after reaching a certain level it becomes stabilised, which can be observed in Fig. 7. The difference between the sumI and sumD controls is subtle, but still sumD results in slightly fewer deaths, as shown in Fig. 7. An interesting observation is that although maxI stabilised the number of symptomatic cases below the peaks in two other strategies, this results in the highest number of deaths.

Fig. 6

NPI strategies aggregated to the country for different objectives for County and Voivodeship models.

Fig. 7

Number of symptomatic cases and deaths for different objectives.

NPI strategies for different regions in the case of the sumI objective are presented in Fig. 8. There are three counties related to three large cities in Poland: Cracow, Wroclaw and Lodz. The least number of initial cases is found in county Cracow, while in Wroclaw 32% more cases are initially recorded. In county Lodz this figure is 56%. It is interesting that the most restrictive NPIs are not applied to the county with the highest ratio of initial cases, but to the county with the least. In such a region the disease is quite suppressed throughout the time horizon, while the most infected region is allowed to reach more than double the number of cases at its peak, although the initial number of cases is not very different in terms of absolute values. We make a note of this interesting observation for future investigation.

Fig. 8

NPI controls and the number of symptomatic cases in different regions for SumI. correspond to NPI level/number of infections at Cracow, Wroclaw and Lodz respectively.

Robustness analysis

Optimisation method

The CMA-ES algorithm we apply to solve the problem is a stochastic search algorithm, so the results obtained may depend on the seed of the pseudorandom number generator used. Due to prolonged computation times, however, it is not possible in practice to analyse multiple model variants along with a large set of different seeds. Fortunately, the optimisation method appears to be robust to different seeds and gives a relatively small spread of solutions. Table 6 presents aggregated computation results for different granularity of decision and various quality criteria for a series of ten algorithm runs with different seeds of a pseudorandom number generator. The deviations of extreme solutions, i.e. smallest and largest, from the averaged solution of ten optimisation runs are given. These deviations do not exceed 0.33%. They reach slightly higher values only for the voivodeship with the sumI criterion, although still significantly below 1%. Therefore, it can be concluded that the optimisation method is resistant to seeds and, due to the long computation times required, a limited number (in most cases just one) of runs of the optimisation process is sufficient.

Table 6

Aggregated optimisation results for ten different seeds: average result and minimum and maximum deviations.

Country
	Min dev.	Average	Max dev.
sumI	0.18%	4 765 694	0.17%
sumD	0.04%	89 026	0.15%
maxI	0.30%	595 351	0.33%
Voivodeship
	Min dev.	Average	Max dev.
sumI	0.88%	4 566 809	0.77%
sumD	0.16%	82 550	0.25%
maxI	0.08%	557 249	0.05%
County
	Min dev.	Average	Max dev.
sumI	0.24%	4 753 168	0.30%
sumD	0.30%	86 534	0.25%
maxI	0.22%	586 253	0.27%

Robustness analysis method

The mobility input data and parameters of the disease development model are estimates only. Therefore, it is important to understand to what extent the NPI policy obtained from optimisation is robust to changes in these data. In order to investigate robustness, we conducted a series of experiments described in the following sections. The single experiment process is illustrated in Fig. 9.

Fig. 9

Schema of robustness analysis method.

Let us denote the estimated input data of the problem introduced in Section 4 as baseline parameters. These include baseline mobility data and baseline disease development parameters. The NPI policy resulting from the solution of the optimisation problem for the baseline data is called the baseline NPI policy. The baseline NPI policy may of course be different depending on the objective assumed.

In a single experiment, we create a scenario by distorting the baseline data, thus obtaining distorted input data. The scenario assumes that the distorted input data is the actual realisation, and the baseline data represents their imperfect estimation. Therefore, the NPI policy obtained as a result of solving the problem for distorted input data is called the proper NPI policy (the policy that should be applied, but its calculation requires perfect knowledge of the input data). In all experiments described in this section we set the maximum number of objective evaluations to 5e5.

Then, the baseline NPI policy and the proper NPI policy are used to simulate the development of the disease under the scenario, thus obtaining the courses of two disease development processes and the values of related objectives. If objectives are equal or close to each other, then applying the baseline policy obtained for imperfect parameter estimation is good enough, and the solution is robust. Taking into account the results obtained in Section 6.1, the solutions can be considered equivalent when the objective does not differ more than the maximum deviation from the mean solution when the optimisation algorithm is run multiple times with a different seed for pseudorandom numbers.

Robustness of mobility data

First, we tested several scenarios for disturbed mobility from the baseline values by a percentage ranging from −100% to 100%. For each scenario, the baseline mobility between each ordered pair of regions has been reduced or increased by a given percentage. Table 7 shows the deterioration of the objective when the baseline NPI policy is applied to distorted mobility, and represents the loss of solution quality. It is noteworthy that the solution is robust to increased mobility in the range up to 50% for the sumI and sumD criteria. The problem with the maxI objective is more sensitive, but up to the level of 20%, the solution is also basically robust. In the case of reduced mobility, the solution remains robust to the level of −40% for the sumI and sumD criteria, and to the level of −10% for the maxI criterion.

Table 7

Deviations of objectives for baseline control with disturbed mobility data.

Mobility disruption	sumI	maxI	sumD
−100%	9.8%	58.8%	6.9%
−99%	7.0%	58.3%	6.1%
−90%	3.7%	19.4%	2.2%
−50%	2.7%	20.2%	1.4%
−40%	0.2%	5.9%	0.2%
−30%	0.2%	3.5%	0.2%
−20%	0.1%	2.0%	0.1%
−10%	0.0%	0.7%	0.0%
＋10%	0.0%	0.3%	0.0%
＋20%	0.0%	0.5%	0.1%
＋30%	0.0%	1.0%	0.1%
＋40%	0.1%	1.8%	0.2%
＋50%	0.9%	7.3%	0.9%
＋100%	1.3%	8.4%	1.3%

We also tested robustness under mobility disturbances introduced randomly. For each ordered pair of regions, , with a baseline number of travellers from to equal to we generate a new number of travellers with normal distribution , where . For such parameters, the solution showed robustness, i.e. the loss of objective quality due to the application of the baseline NPI policy in relation to the optimised NPI policy for disturbed mobility was below 0.5%, which is a negligible value. This is probably related to the fact that the exact values of the number of travellers are not so important as their relations. In order to illustrate that, let us consider a big city with a large stream of travellers to another big city and a small stream to many small areas. Significant distortion will still keep large streams of travellers between large cities and small streams to other areas. Thus, random deviations around the baseline mobility values are not significant for the obtained result, as opposed to a general underestimation or overestimation.

Robustness of disease spread model parameters

There is a number of parameters in the disease spread model that can be hard to estimate. Similar to our research in the case of mobility data, we consider some disturbances in the parameters, and we compare the objective optimised for disturbed parameters with the objective obtained from the disease spread simulation under disturbed parameters for the baseline NPI policy. We consider each parameter separately, except for gammas, which are perturbed together.

Table 8 shows the impact of perturbations of contact transmission rates and . Optimised NPI policy weakly depends on contact transmission rates and for objectives sumI and sumD. The only significant impact is visible for lower than the baseline value. For the objective , the impact of perturbations of contact transmission rates is significant and slowly increases with increasing disturbances. Although the solution can be improved by several percent, the baseline NPI policy still remains quite good for the disturbed parameters.

Table 8

Deviations of individual criteria for control with disturbed .

β	sumI	maxI	sumD
0.1	0.4%	2.3%	0.4%
0.15	0.3%	15.7%	0.3%
0.2	Baseline value
0.25	0.0%	0.7%	0.0%
0.3	0.0%	2.0%	0.0%
0.35	0.0%	9.5%	0.1%
0.4	0.1%	3.2%	0.1%
0.45	0.1%	3.8%	0.2%
0.5	0.2%	4.1%	0.2%
μ	sumI	maxI	sumD
0.1	0.9%	5.7%	0.8%
0.15	4.3%	7.7%	2.9%
0.2	Baseline value
0.25	0.1%	3.0%	0.1%
0.3	0.3%	4.8%	0.3%
0.35	0.4%	5.6%	0.4%
0.4	0.4%	6.4%	0.3%
0.45	0.4%	6.3%	0.3%
0.5	0.4%	6.7%	0.3%

The model seems to be the most sensitive to perturbations in recovery rates and . We assumed the following four cases of recovery times , where a pair means days to full recovery for asymptomatic cases, and days to full recovery for case of symptomatic ones. The results are presented in Table 9. The greatest sensitivity of the model is towards the reduction of recovery time, and in general, the maxI model is less robust.

Table 9

Deviations of individual criteria for control with disturbed and .

γ_A	γ_I	sumI	maxI	sumD
0.14 (1/7 days)	0.09 (1/11 days)	4.1%	9.2%	3.4%
0.10 (1/10 days)	0.07 (1/15 days)	Baseline value
0.07 (1/14 days)	0.05 (1/21 days)	0.9%	2.9%	1.0%
0.05 (1/21 days)	0.03 (1/30 days)	0.0%	2.2%	0.1%

While the baseline incubation time (equal to ) was originally assumed to be 5 days, we tested the model for 3, 4, 6, and 7 days. In the case of the sumI and sumD criteria, no significant impact on optimal NPI policy has been observed. The results for maxI are presented in Table 10. While 1 day is safe, bigger disturbances bring loss of objective quality.

Table 10

Loss of objective quality for different .

σ	maxI
0.33 (1/3 days)	2.4%
0.25 (1/4 days)	0.9%
0.20 (1/5 days)	Baseline value
0.17 (1/6 days)	0.8%
0.14 (1/8 days)	1.6%

Computations for different values of the fraction of individuals who become symptomatic show that the baseline NPI policy is adequate for disturbed . The loss of quality of the objective is no greater than 0.5% in each case. This is reasonable since, in our model, symptomatic and asymptomatic cases are equally contagious. Therefore, some significant difference could appear when the change in would be combined with changes in other disease-spreading parameters.

The size of the model parameter space is too large for us to research the effect of any possible combination of parameter disturbances in practice. Single-parameter studies show that our approach seems to be quite robust to parameter changes. We suggest continuing the research if evidence points to any possible parameter inaccuracy scenarios of practical importance.

Conclusions

We introduced an NPI optimisation model based on a network model of disease spread. Our model takes into account the transmission of a disease by individuals travelling to other regions, including transmission among visitors coming to the same region from different origins, which is neglected in many papers. The resulting model is a non-linear, potentially large-scale, programme that is hard to solve. Therefore, we used evolutionary strategy, and we showed that it is possible to decrease the number of symptomatic cases, deaths and daily peaks of cases compared to simple heuristics. Depending on the objective, the obtained control differs significantly, so the objectives are partially contradicted.

Our results show that, in general, it is better to use a more detailed network model. The best results for maxI and sumD were obtained for the County model. When a more, however, detailed model is used the complexity of optimisation tasks is burden some, as can be seen in the case of sumI, for which the County model was not the best, and the optimisation process was stopped after the assumed time of computation. Detailed analysis showed that the objective function was still decreasing at that stage, which meant that it could still be improved.

We assess that with the optimisation technique that we used, a good trade-off between optimisation capabilities and the quality of the spread model lies somewhere between the Voivodeship and County models, which means between 16 and 380 nodes. The Community model was too big to be solved with the evolutionary heuristic we used, but was still promising in terms of the quality of the simulation. Therefore, there is strong need for more efficient optimisation techniques for spatial SIR-like models that could allow for the solution of large-scale problems.

NPIs have wide-ranging effects on the health, economics and social well-being of populations, which has resulted in growing pandemic fatigue and a drop in adherence to NPIs compared to when they were initially implemented. The overall effectiveness of NPIs shows a downward trend between waves, possibly due to policy fatigue (Ge et al., 2022). This is also due, among other things, to the lack of information supported by scientific research on how the pandemic will develop in the long term and what NPIs policy-makers are planning to introduce to overcome it. Future uncertainties resulting from NPIs, such as when and how the measures would be eased and the impact on job prospects, was reported as a major stressor on mental health (Schneiders et al., 2022). Uncertainties about the future resulting from NPIs also strongly affect the ability to make business decisions. Our model allows the determination of how the NPI policy regarding COVID-19 should be shaped in the longer term. Hence, our model is a tool that supports policy-makers so that they can inform the public about predictions regarding long-term NPI policy plans.

CRediT authorship contribution statement

Mariusz Kaleta: Conceptualization, Methodology, Software, Investigation, Data curation, Writing – original draft. Małgorzata Kęsik-Brodacka: Resources, Writing – reviewing and editing. Karolina Nowak: Resources, Writing – review & editing. Robert Olszewski: Visualization, Data curation, Writing – review & editing, Funding acquisition. Tomasz Śliwiński: Conceptualization, Methodology, Software, Writing – reviewing and editing. Izabela Żółtowska: Conceptualization, Methodology, Writing – reviewing and editing.