Abstract
Many countries are currently dealing with the COVID-19 epidemic and are searching for an exit strategy such that life in society can return to normal. To support this search, computational models are used to predict the spread of the virus and to assess the efficacy of policy measures before actual implementation. The model output has to be interpreted carefully though, as computational models are subject to uncertainties. These can stem from, e.g., limited knowledge about input parameters values or from the intrinsic stochastic nature of some computational models. They lead to uncertainties in the model predictions, raising the question what distribution of values the model produces for key indicators of the severity of the epidemic. Here we show how to tackle this question using techniques for uncertainty quantification and sensitivity analysis.
We assess the uncertainties and sensitivities of four exit strategies implemented in an agent-based transmission model with geographical stratification. The exit strategies are termed Flattening the Curve, Contact Tracing, Intermittent Lockdown and Phased Opening. We consider two key indicators of the ability of exit strategies to avoid catastrophic health care overload: the maximum number of prevalent cases in intensive care (IC), and the total number of IC patient-days in excess of IC bed capacity. Our results show that uncertainties not directly related to the exit strategies are secondary, although they should still be considered in comprehensive analysis intended to inform policy makers. The sensitivity analysis discloses the crucial role of the intervention uptake by the population and of the capability to trace infected individuals. Finally, we explore the existence of a safe operating space. For Intermittent Lockdown we find only a small region in the model parameter space where the key indicators of the model stay within safe bounds, whereas this region is larger for the other exit strategies.
Author summary Many countries are currently dealing with the COVID-19 epidemic and are looking for an exit strategy such that life in society can return to normal. For that purpose computational models are used to predict the spread of the virus and to assess the efficacy of policy measures before putting them into practice. These models are subject to uncertainties (due to, for instance, limited knowledge of the parameter values), which can lead to a large variability in model predictions. It is therefore fundamental to assess which range of values a model produces for key indicators of the severity of the epidemic.
We present here the results of the uncertainty and sensitivity analysis of four exit strategies simulated with an individual-based model of the COVID-19 transmission. As key indicators of the severity of the pandemic we consider the maximum number of cases in intensive care and the total number of intensive care patient-days in excess. Our results show the crucial role of the intervention uptake by the population, of the reduction in the level of transmission by intervention and of the capability to trace infected individuals.
Introduction
Many countries are currently dealing with the COVID-19 epidemic and are searching for an exit strategy such that life in society can return to normal. However, in absence of an effective curative treatment and, until recently, of an effective vaccine, non-pharmaceutical interventions have been used to keep case numbers as low as possible. In the past there have been numerous other epidemics during which government actions were required to protect the population. Examples are the influenza pandemic (also known as Spanish flu) in 1918 or the more recent Mexican flu (or swine flu) in 2009 [1]. When an infectious disease outbreak occurs, governments rely on computational models to predict the spread of the disease and to explore the potential impact of interventions [2–4]. Thus, computational models enable decision makers in government and public health institutions to assess the efficacy of policy measures before actual implementation.
Since computational modeling of the epidemic has come to play a significant role for informing policy, it is important to assess the uncertainties of the models and of their predictions. Such uncertainties can stem, for instance, from limited knowledge about the values of the input parameters or from the intrinsic stochastic nature of part of the computational models. The presence of any of these uncertainties leads to uncertainties in the model predictions. A central question is therefore what distribution of values is produced by the model for key indicators of the severity of the epidemic and of the intensity of interventions.
In this study we present results from an analysis of uncertainties and sensitivities of an agent-based model of the COVID-19 epidemic [5,6], obtained using techniques of Uncertainty Quantification (UQ) and Sensitivity Analysis (SA) [7,8]. UQ is an area of mathematics dealing with propagation of uncertainties from model input to model output, and establishing model input uncertainties from observation data. Let us denote with ℳ the computational model and let X and Y be its input parameters and output quantities, respectively, such that Y =ℳ (X). X has probability distribution p (X), reflecting its uncertainty. A central aim of UQ is to determine the probability distribution of Y, i.e. p (Y), given p (X). With SA, one determines what fraction of the uncertainty of Y (e.g. fraction of the total variance of Y) is due to individual elements of the vector X, given p (X). Thus, one searches which input parameters generate the most (or the least) variance in Y. It is frequently the case that a large part of the variance of an output quantity is due to a subset of the input parameters involved.
The UQ and SA frameworks employ a probability distribution - as opposed to a single value - to describe each input parameter. These distributions are typically determined from available data or from expert knowledge. With proper probability distributions assigned to each input parameter, the UQ and SA frameworks can be used to assess whether exit strategies are robust. A useful criterion for robustness is that, given the model input distributions, the 95th percentile of a chosen output quantity - usually termed Quantity of Interest (QoI) - remains below a critical threshold. Examples of QoI in the context of epidemic modeling are the peak number of COVID-19 patients in intensive care (IC) units and the total number of fatalities due to COVID-19. For other QoIs, e.g. life years gained, a more relevant criterion is that the 5th percentile stays above a minimum value. This can be dealt with in an analogous manner.
The aim of this work is to perform a model-based quantitative analysis of uncertainties and sensitivities for the computational representations of four exit strategies, and assess the uncertainties in model simulation results. We show how such analysis can be performed by means of computational methods and concepts from the fields of uncertainty quantification and sensitivity analysis, and what kind of insights can be obtained. In order to perform our analysis we consider the spread of the COVID-19 disease in the Netherlands in the context of an open-source agent-based model with geographical stratification [5,6]. We note that the framework of our analysis is not restricted to this specific model, but can also be applied to other strategies, to more complex models and to other epidemics than COVID-19.
Methods and model
In what follows, we provide an overview of the computational model used in our analysis and the exit strategies implemented in the model. Next, we summarize some key concepts of uncertainty quantification and describe different sources of uncertainties. We discuss also the chosen SA method and the quantities of interest selected for our study. We conclude this section with a description of the computational UQ and SA framework that we employ.
Computational model
We employ the publicly available virsim model [6], which is a stylized representation of the COVID-19 epidemic with geographical stratification of both transmission and interventions (i.e., a meta-population model). More specifically, it is an agent-based model with a geographical structure defined by means of clusters (representing for instance towns and villages) and superclusters (e.g. provinces or other administrative units for which policy decisions are made). Individuals are separated into susceptible (S), exposed (E), infectious (I) and recovered (R), and life-long immunity is assumed. Heterogeneity is introduced by varying individual contact rates. For a more detailed model description we refer the reader to de Vlas & Coffeng [5]. For this study, we employed an updated quantification of the model [9]; see Table 1 in the main text and S1 Table for an overview of the parameter values.
We consider four exit strategies: Flattening the Curve (FC), Contact Tracing (CT), Intermittent Lockdown (IL) and Phased Opening (PO). Each strategy is part of the model implementation (see S1 Appendix for the computational details) and is defined by a unique set of parameters. More details about the strategies parameters are given in a later section, where we discuss also the input distributions that have been considered to account for the respective uncertainties.
Since this study is of conceptual nature and does not aim to model real-world scenarios as accurately as possible, we simulate a population of 1 million individuals and we focus on the first year after the implementation of a strategy (the idea being that governments will regularly evaluate and adapt the strategy in use). All simulated intervention strategies are preceded by a period of lockdown, analogous to the situation in the Netherlands from March to May 2020. During this period, we assume that an intervention package is implemented that heavily reduces transmission. The model parameter for this effect (intervention_effect) is expressed in terms of the level of transmission that the intervention package still allows, relative to a situation without any interventions (i.e., 100% minus the reduction in the overall transmission rate).
Based on a previous quantification of the model using extensive data on the Dutch COVID-19 outbreak [9], we assume that the lockdown reduces transmission to 30% on average for one week and then to 15% on average after the Dutch government introduced stricter measures. In accordance with the Dutch situation, the effect of the lockdown is simulated for 60 days, followed by a period of partial lockdown (assumption: intervention_effect = 25%), analogous to the re-opening of schools in the Netherlands [9]. After another 30 days, one of the exit strategies takes effect in the simulation. Overall we simulate a period of 550 days corresponding to roughly 1.5 years; see S1 Appendix for details about how this is implemented in the model.
Uncertainty and its quantification
Uncertainties in model output can arise from different sources; we discuss here four main types. The first is parameter uncertainty, referring to uncertainties in model parameters whose values can be set directly by the model user via the inputs of the computational model. An example is the reduction of the transmission rate due to the introduction of an intervention. For the COVID-19 epidemic, a significant decrease in the transmission rate following the adoption of such measures has been recorded [10–13], however the magnitude of this reduction is uncertain.
The second type of uncertainty, which we call intrinsic uncertainty, arises when a computational model is inherently stochastic. In epidemiology, many models are agent-based and possess inherent stochasticity, for instance in the randomized interactions between agents. Model users often have little control over such inherent stochasticity as they can typically only set the seed of the random number generator at the start of a simulation [14].
The third type is model-form uncertainty, referring to uncertainty or errors in the structure of the model itself (e.g. due to transmission mechanisms not represented in the model). This type of uncertainty cannot be analyzed by changes in the model inputs but requires a comparison either with independent observation data or with other models - as done for example in climate science [15].
Lastly, initial condition uncertainty is due to the inaccuracies in the specification of the initial state of the model (i.e., the state of the simulated population at the start of the model run). Since we consider here model outputs which are independent from the specific timing of e.g. epidemic peaks, this type of uncertainty is not important in our analysis.
In this study we analyze the parameter and intrinsic uncertainties by means of non-intrusive UQ methods, which means that we treat the computational model (i.e., the virsim model) as a black box. We extend the notation from the introduction to Y = ℳ (X, r), where r denotes the seed of the random number generator. By random sampling of the parameter vector X from its distribution p (X) and altering the seed r, followed by executing the model ℳ, we create random samples from the probability distribution of Y. Thus, we can probe the uncertainty of Y and estimate its distribution by repeated model executions.
Sensitivity analysis
In order to assess which parameters create the most uncertainty in the model output Y under the different strategies, we compute the Sobol indices [16]. This is a form of global sensitivity analysis focusing on the variance of Y. Assuming mutual independence among the input parameters, Var(Y) is decomposed into fractions which can be attributed either to a single input parameter (first order Sobol indices) or to a set of parameters (higher order Sobol indices). Given p (X), the first order Sobol index for the i-th input parameter (Xi) is defined as Si := Vari[𝔼∼i(Y |Xi)] / Var(Y), as explained in S2 Appendix. We use the notation 𝔼i for expectation over Xi and 𝔼∼i for expectation over all X1, X2, … except Xi (and similarly for Vari and Var∼i). If Si is close to 1, it means that the variance of Y is almost entirely due to the variance of Xi. The overall effect on the model output of all parameter combinations involving Xi is given by the total Sobol index .
Regarding the assumed mutual independence of the input parameters, there might be dependencies when parameters are actually estimated from data. However such dependencies are not implemented in the computational model, hence the selected input parameters (see below for more details) are effectively treated as mutually independent. We refer to S2 Appendix for more details about the theory of the Sobol indices.
Quantities of Interest
In our analysis we consider the model predictions for the number of incident and prevalent individuals in the population that require IC admission. As IC capacity is limited, the question whether the IC capacity will be exceeded (and if so, by how much) according to model predictions is clearly important. To investigate this, for each simulation we consider two quantities of interest (QoI):
the maximum of the moving average of the prevalent cases in intensive care (averaging window = 30 days)
the total number of IC patient-days in excess of IC bed capacity (referred to as “IC patient-days in excess” from hereon).
The first QoI shows the peak value of the number of COVID-19 patients in IC units, giving an indication of the intensity of an outbreak. We apply a moving average to focus on longer-term trends, filtering out short-term “noisy” variations. When analyzing robustness, a natural threshold for this QoI is the available IC capacity, which may vary from country to country and from month to month. In our analysis we assume the IC capacity to be constant in time and we consider the Netherlands as reference country. De Vlas and Coffeng [5] calculate that during the first epidemic wave in 2020 the Dutch maximum IC capacity for COVID-19 patients was 109 IC beds per million inhabitants.
The second QoI quantifies by how much the total IC capacity is overburdened. It is defined as where p(t) is the number of IC prevalent cases at day t∈ {0, 1, 2, …, T}, p* is the maximum IC capacity and the summation runs from the start time of the simulation (t = 0 day) to the final time (t = T days). It shows the extent to which the government may have to rely on other countries, e.g. Germany, to take in IC patients from the Netherlands. Setting a threshold for this QoI is complicated, as many economical and political factors come into play here.
We note that the two QoIs as defined here do not depend on time. The first QoI is defined as a maximum over the time interval of simulation, whereas the second QoI is defined as a summation over the same time interval. Thus, a single model simulation yields a single, time-independent value for each QoI.
Policy-related uncertainties
The implementation of a strategy in the computational model is determined by a set of input parameters for the model. Some of the parameters can be controlled (up to a certain degree) by policy makers. These are parameters related to social aspects of the population - e.g. social distancing - and to the availability of resources, e.g. to test and track infected individuals. They can be controlled to some extent by government authorities through imposing stricter (or less strict) rules for e.g. social distancing, or through making more (or less) resources available. We call them policy-related parameters and we refer to the uncertainties generated by these parameters as policy-related uncertainties.
Below we provide details about which policy-related parameters are treated as uncertain, together with the rationale behind the probability distributions chosen for them. The distributions, together with their mean and 95%-confidence intervals are provided in Table 1; their plot instead is provided in S1 Fig.
We use Beta distributions for those parameters that are naturally restricted to values on a finite interval, such as those representing probabilities or percentages. We remark that this allows us to use distributions that are very different from the uniform distribution. With the hyperparameters used here, the probability density of the Beta distribution decreases to zero towards the boundaries of the support. By contrast, the uniform distribution attributes equal probability to any value within the support of the distribution. For parameters defined as a time scale or a rate we choose the Γ distribution as it has one unique peak (with the hyperparameters used here) and a semi-infinite domain (all non-negative values).
Characteristics of the four exit strategies
Flattening the Curve
The main idea of FC is to gradually resume activities that have been interrupted by the lockdown such that the virus spreads among the population at a lower pace while keeping the pressure on the healthcare system manageable. Here we model only the first intervention after the lockdown with the idea that the evolution of the pandemic can be appropriately reconstructed up to (almost) present time, such that the focus is on the next step to take. In the computational model, this strategy is governed by two parameters:
the effect of the intervention (model parameter intervention_effect). The effect of the intervention package is expressed in terms of the level of transmission it still allows, relative to a situation without any intervention (i.e., one minus the reduction in the overall transmission rate). This parameter is supposed to allow a slight increase in the average contact rate with respect to the situation of general lockdown with schools open. We sample it from a Beta distribution;
the uptake of the intervention by the population, model parameter uptake, which we draw from a Beta distribution. We assume that the majority of the population would adhere to the introduced measures, hence the bulk of the distribution is moved towards values close to 1.
Contact Tracing
This strategy consists of tracking potentially infectious contacts such that only infected people drastically reduce their interactions with others, thus limiting the spreading of the virus while allowing events and activities to take place nonetheless. Three factors characterize such an approach:
the delay between becoming infectious and being identified as such (if at all); represented by the inverse of the model parameter trace_rate_I. A fast identification is technologically very challenging, hence we give higher probabilities to identifications requiring more than 20 hours. We employ a Γ distribution for the inverse of the delay;
the probability of an infected contact to be identified before the person turns infectious; represented by the model parameter trace_prob_E. Such a classification is challenging as tracking data might be incomplete or there might not be enough capacity to process them. Therefore we choose a Beta distribution with relatively large variance and mean shifted towards lower values;
the quality of the isolation and its effect on transmission; represented by the model parameter trace_contact_reduction. We assume that whoever is being tracked would adopt every possible measure to avoid transmitting the virus further. For this parameter we use a Beta distribution.
Intermittent Lockdown
The idea here is to alternate periods of lockdown and periods of complete opening at preset intervals (at the time of writing the virsim model did not include a dynamic trigger for the start of the lockdown/opening based on, e.g., the number of infected). When such a strategy is to be modeled, the following parameters have to be chosen:
the duration of the lockdown and of the following lift of measures represented by the model parameters lock_length and lift_length. Since a long lockdown might have a heavy impact on the psychological health of the population with more people, for instance, becoming depressed, we try to keep a balance between the lockdown and the lift period. We consider Γ distributions both for the lockdown and for the lift periods;
the effect of the lockdown (model parameter lock_effect) expressed in terms of the level of transmission it still allows, relative to a situation without any intervention. We assume it to have (on average) an effect in contact reduction similar to a general lockdown with schools open. We sample it from a Beta distribution;
the uptake of the intervention by the population. For this parameter we use the same Beta distribution assumed in case of FC.
Phased Opening
The approach of the Phased Opening strategy is to release the general lockdown at the regional level and distribute the patients in need of hospital care at the national level. In this way the spread of the virus is limited to a smaller portion of the population and the burden on the healthcare system is reduced [5]. In our simulations we consider a number of phases equal to the number of superclusters, so at each phase only one supercluster lifts the lockdown. Other important features defining the strategy are:
the time interval between one phase and the next, determined by intervention_lift_interval. A too short period would lead the strategy to be rather inefficient and to a strong overburden on the healthcare system in a short time. On the other hand a too long period would lead to a rather long waiting for some regions and hence to growing discontent among the population. We try to find a compromise in the choice of the Γ distribution for this parameter;
the effect of the measures in the areas where the lift has not been applied yet, represented by the model parameter pl_intervention_effect_hi and expressed in terms of the level of transmission it still allows, relative to a situation without any intervention. In these areas we assume that the reduction in the average contact rate is analogous to the reduction obtained during the periods of lockdown in case of the IL strategy, hence we use the same Beta distribution;
the uptake of the intervention by the population, model parameter uptake, which is sampled from a Beta distribution, as in the cases of FC and IL.
Other uncertainties
Besides the parameters that can be influenced by policy, discussed in the previous section, the computational model has other parameters that we include in our analysis. Examples are the reproduction number of the virus and the duration of infectiousness. Because these parameters have to be estimated from data, their values can be rather uncertain, especially when available data is scarce (for instance in the early stages of an epidemic with a new virus). For some parameters, care has to be taken when sampling them from probability distributions, since certain basic characteristics, like the doubling time, have to match the time evolution of the epidemic as captured in the (real-world) observation data.
Altogether, given the setup of the virsim model, we consider the following non-policy-related uncertain parameters:
the average duration of infectiousness (model parameter avg_duration_infectiousness) has been estimated to be about 5 days. To account for errors and a conceivable lack of data in the estimation process, we sample it from a Γ distribution;
the reproduction number R0 is determined in the model as the product between the average duration of infectiousness and the contact rate. By default it is set at R0 = 2.5 and in what follows we account for uncertainties on how the virus spreads among the Dutch society. Thus we draw R0 from a Γ distribution and derive the average contact rate as the ratio between R0 and the average duration of infectiousness;
intervention_effect_var represents the inter-individual variation in the effect of the intervention in case of uptake (in the default setup no inter-individual variation is allowed). Such variations might be due, for instance, to the presence of children in the household. intervention_effect_var can assume any positive value and, for simplicity of representation, it is modeled as 1/x where x is sampled from a Γ distribution. Please note that even though the 95%-confidence interval for x itself is narrow, the resulting interval for 1/x is quite large;
in order to consider uncertainties in how the incubation period varies among the population, we sample the shape parameter of the distribution of exposed_time from a Γ distribution.
The virsim model is stochastic, and as discussed in an earlier section, this intrinsic stochasticity is a source of uncertainty of the model output. It is not possible to choose the probability distribution of the internal (or latent/hidden) random variables of the model, or to set them by hand to specific values according to some sampling plan. The only form of control as model user is to pick the random seed for the pseudo-random number generator at the start of the simulation. We draw the seed from a discrete uniform distribution with support between 16384 and 65536. This amounts to independent random sampling of the internal random variables of the virsim model.
Furthermore, when using models with a geographical structure like the virsim model, uncertainties can arise from the level of geographical mixing that is allowed in the model. For sake of simplicity we do not consider such uncertainties in the present study, but they should be considered in more comprehensive studies if these are intended to inform policy makers.
We report in Table 1 the assumed input distributions, their mean and 95%-confidence intervals reflecting parameter uncertainty for model parameters that were considered in the UQ analysis. The plot of the distributions themselves is provided in S1 Fig. Finally, in S1 Table we list the model parameters that have fixed values (i.e., the parameters not considered as uncertain) in this study.
UQ and SA computational framework
If the dependence of the QoIs on the parameters is smooth and the number of uncertain parameters is not too high, the propagation of uncertainties from parameters to QoIs can be assessed with techniques such as Polynomial Chaos Expansion and Stochastic Collocation [7, 8, 17], that allow to obtain the relevant information with relatively little computational time. However, due to the presence of intrinsic uncertainty, the continuous relation between QoIs and inputs is lost as the QoIs do not depend continuously on the random seed. When this type of uncertainties is present, the aforementioned methods like Polynomial Chaos Expansion are not suitable. Therefore, we resort to Monte Carlo (MC) sampling as it is guaranteed to return reliable results (albeit with a low convergence rate). We base our UQ results on the evaluation of an MC ensemble with 1000 simulations. More advanced techniques meant to overcome the issue of the internal stochasticity are currently being developed, see for instance [18] but require further investigation before application to high-dimensional systems.
For the SA, we compute the first order Sobol indices with the cost effective algorithm of Saltelli [19] (see S2 Appendix), and their 95%-confidence intervals are computed via bootstrapping. By performing an intelligent sampling of the input parameters, this algorithm reduces the computational cost of the calculation of the first order Sobol indices from M 2p to M (p + 2) model simulations, where M is the number of MC samples and p is the number of uncertain parameters (including the random seed). For more details about the practical implementation of the Saltelli algorithm please see [8]. Due to the specific sampling required by the Saltelli algorithm, the SA cannot be performed on the same set of data employed for UQ. Therefore we run a separate set of simulations, and use M = 2000 MC samples per parameter involved in the implementation of the strategy (SA is in general computationally more demanding than UQ). This results, e.g., in case of FC in a total of M (p + 2) = 2000(3 + 2) = 10000 simulations, where we considered only the random seed and the policy-related parameters.
Sampling and post-processing analysis are done using the Monte Carlo sampler of the publicly available Python library EasyVVUQ [20–22]. We run the required model simulations in parallel on a supercomputer at the Poznan Supercomputing and Networking Center. For the job submission to the supercomputer we use the FabSim3 [23,24] and QCG-PilotJob [25] packages. FabSim3, QCG-PilotJob and EasyVVUQ are all part of the open source VECMA Toolkit (VECMAtk) [26, 27]. The codes employed for this study can be found on GitHub [28].
Results
To provide intuition for the behaviour of the model and the QoIs defined earlier, we show in Figure 1 time series of the number of prevalent cases in IC and of the number of IC patient-days in excess, for four simulations of the virsim model under the Phased Opening strategy. For each simulation, different policy-related parameters are used. It can be seen in the left panel of Fig. 1 that the number of prevalent cases in IC has noisy perturbations on top of the main signal. As discussed earlier, we apply a moving average (averaging window = 30 days) to filter out these short-term noisy perturbations. The resulting smoothed time series are shown in the central panel of Fig. 1. The peak value reached in each smoothed time series is our first QoI. In the right panel of Fig. 1 we show the cumulative number of IC patient-days in excess as a function of time.The value reached at the final time, t = T, is our second QoI.
Probability distributions of the QoIs
As first step of our UQ analysis, we construct the empirical cumulative distribution functions (cdfs) of our QoIs. For any threshold value q*, the cdf gives the probability that the QoI remains below (or at) that threshold, i.e. it gives P(QoI ≤ q*) for the adopted strategy (and given the specific input parameter ditributions). In Fig. 2 we report the resulting cdfs for the four selected strategies both with and without uncertainties in non-policy-related parameters (in the latter case, we fix these parameters at the mean values listed in Table 1, except for shape_exposed_time which is set to 20, see also S1 Table, and for the seed which is never fixed). We further display the 95%-confidence interval of the empirical cdfs using the Dvoretzky-Kiefer-Wolfowitz inequality [29,30], which gives an indication on how reliable the empirical cdfs are.
We observe the important fact that, with the given distributions of the input parameters, none of the analyzed strategies is robust. The probability that the number of prevalent patients in intensive care is larger than the IC capacity is rather high and only Contact Tracing gets close to a probability of 50%. This shows that the assumed input distributions for Flattening the Curve, Intermittent Lockdown and Phased Opening correspond to interventions that are not sufficiently restrictive to stay below the threshold, as far as the model is concerned.
We note that the cdfs of the first QoI for CT and IL increase more gradually compared to FC and PO (implying that the variance of the first QoI is larger for CT and IL than it is for FC and PO). Furthermore, it can be seen that the shape of the cdfs is only weakly affected by non-policy-related parameters. Therefore, when searching for the parameters responsible for most of the output variability, i.e. for the sensitivity analysis, these parameters might be kept fixed to reduce the computational burden. They should be included however when determining the minimum level of e.g. intervention or uptake, required for the QoIs to stay below their threshold with 95% probability.
Sensitivity analysis
In Fig. 3 we report, for the first QoI, the 95%-confidence interval of the first order Sobol indices for the policy-related parameters and for the seed. This information allows us to identify which factors of the strategies most affect the model output, such that measures targeting these factors can be adopted in order to make the strategy robust. The Sobol indices for the second QoI show qualitatively similar results and are therefore not reported here.
From the width of the intervals, it can be seen that the number of MC samples (M = 2000) is too small to estimate the Sobol indices very accurately. Nonetheless a qualitative ranking of the parameters with respect to the generated output variance can be made, based on the results in Fig. 3. It is therefore possible to distinguish intervention_effect, trace_rate_I, uptake and pl_intervention_effect_hi as the parameters responsible for most of the output uncertainty, and should thus be targeted (e.g. through policy measures) in order to obtain more desirable outcomes. For more accurate estimates of the Sobol indices, M must be increased substantially, thereby greatly increasing also the computational cost. We refrain from doing so in this study.
The uptake by the population of the interventions plays a crucial role whenever this parameter is part of the strategy. In these strategies (i.e. FC, IL, PO), the uptake parameter is responsible for at least 30% of the QoI variance and, in case of Intermittent Lockdown it is responsible for about 70% of the variance. The intervention_effect, lock_effect and pl_intervention_effect_hi parameters are also important, although they are not always the main driving factor. For Flattening the Curve its Sobol index is higher than that of uptake, whereas for Intermittent Lockdown uptake has a much higher Sobol index (indicating that is it more important) than lock_effect. The population uptake and the effect of the lockdown seem to be of similar importance for the Phased Opening strategy.
In case of Contact Tracing circa 50% of the QoI variance is determined by the delay between becoming infectious and being identified as such (if at all), i.e. the inverse of the rate per day at which infected individuals are being traced (parameter trace_rate_I). After infected agents have been identified, the reduction of their average contact rate is also important.
The probability of tracing exposed individuals (CT strategy) and the interval between subsequent lifts (PO strategy) give only a small contribution to the model output variability. The intrinsic stochasticity of the model instead does not induce much variability in the model output as the 95%-CI of its Sobol index always includes 0 and does not take values above 5%. Similar conclusions hold for the lengths of the lockdowns and subsequent lift periods in the IL strategy.
The total Sobol indices give qualitatively the same outcomes but highlight higher order interactions involving lock_effect and trace_prob_E (see S2 Fig). The differences between the first order and the total Sobol indices suggest that the most important second-order interactions appear to be those between lock_effect and uptake in the IL strategy, and between trace_prob_E and trace_rate_I in CT. However, to assess these interaction effects in detail, more MC samples are needed.
Safe operating space
Given the knowledge on the main driving factors of each strategy, we want to know which combinations of values for these input parameters result in, for instance, a number of prevalent cases in IC (first QoI) below the IC capacity. This information can be used to devise policy measures that would effectively move the corresponding input distributions towards an area in the parameter space where the strategy is robust, i.e. towards the safe operating space. The information can be obtained by means of a scatter plot of the QoIs as functions of the values of the two or three input parameters with the highest Sobol indices.
In case of Flattening the Curve we visualize the two selected QoIs as functions of the most important parameters according to the Sobol indices, i.e. intervention_effect and uptake, see Fig. 4. The first QoI stays below the IC capacity for high population uptake and strong intervention (recall that stronger intervention is associated with lower values of intervention_effect). If uptake is less than circa 85% or if intervention_effect is higher than 35%, an exceedance of the IC capacity is to be expected in the model.
The two main drivers of the strategy Contact Tracing are trace_rate_I and trace_contact_reduction, while trace_prob_E has higher relevance when considered in combination with the other parameters. Therefore, we ignore the seed parameter as the variance it induces is small (the 95%-CI of its Sobol index being very tight around 0) and split our data into quartiles (very low, low, high and very high) of trace_prob_E. The scatter plots of the resulting sub-sets of data of the first QoI are provided in Fig. 5. The same analysis for the second QoI reveals the same qualitative results and is therefore not reported here.
At low values of trace_rate_I (i.e. in case of a long delay in the identification process), the first QoI is predicted to be above the threshold independently from trace_contact_reduction. This is slightly mitigated in the plots corresponding to the high and very high quartiles of trace_prob_E, revealing the important interaction between these two parameters captured by the total Sobol index. It is therefore paramount to have an efficient tracing system for this strategy to be effective.
Out of the five considered parameters of the strategy Intermittent Lockdown, only two are important according to the estimated Sobol indices. These are uptake and lock_effect. The scatter plots of the QoIs as functions of these two parameters (see Fig. 6) confirm that high levels of uptake are necessary when adopting this strategy. They also show that, as the level of transmission allowed by the lockdown increases (i.e. higher values of lock_effect), high levels of uptake are not enough to keep the first QoI below its threshold. Very little variability is induced by the other parameters.
Similar to Contact Tracing, the Phased Opening strategy has two parameters responsible for most of the output variance: uptake and pl_intervention_effect_hi. The interval between consecutive phases, i.e. intervention_lift_interval, has a small Sobol index, whereas the Sobol index of the random seed is very close to zero. We therefore apply a similar approach: we ignore the random seed and divide the data into quartiles of intervention_lift_interval. We show the scatter plots of the resulting sub-sets of the first QoI in Fig. 7. The same analysis for the other QoI reveals the same qualitative insights and is therefore not reported here.
There is a clear gradient visible in Fig. 7 with lowest QoI values obtained for high values of uptake and for low values of pl_intervention_effect_hi (corresponding to more restrictive lockdowns). Higher QoI values are obtained when either the population uptake or the rigor of the lockdown decreases. These two factors alone though are not enough to ensure that the IC capacity is not exceeded. In fact all of our simulations in the quartile corresponding to very low values of intervention_lift_interval go beyond the IC capacity (top left panel in Fig. 7). By increasing intervention_lift_interval, more simulations stay below the IC capacity threshold.
Discussion
The aim of this work was to perform a model-based quantitative analysis of the uncertainties and sensitivities of a number of selected exit strategies for the COVID-19 epidemic. We discussed how such an analysis can be approached using methods and concepts from the field of Uncertainty Quantification and Sensitivity Analysis. We demonstrated computational techniques that can be employed to assess the uncertainties and identify the sensitivities of each strategy. In particular, we examined the empirical cumulative distribution function of two quantities of interest obtained from the model output: the maximum number of prevalent cases in IC and the total amount of IC patient-days in excess of IC bed capacity. We also identified the input parameters responsible for most of the output uncertainty (namely intervention_effect, trace_rate_I, uptake and pl_intervention_effect_hi) via the computation of the Sobol indices of variance. Lastly we explored the shape of the safe operating space in the parameter space of the different strategies by means of scatter plots.
Given the probability distributions that we chose for the uncertain model parameters, the Contact Tracing strategy appears to be the most effective in the model (see Fig. 2), but it requires high tracing capacity. Furthermore, the output uncertainty is highly driven by the uncertainty in the population uptake of, and adherence to interventions (see Fig. 3), with the number of IC patients increasing as the value of uptake decreases and vice versa. The effect of interventions and the tracing capability also play a crucial role, whereas the intrinsic stochasticity of the model always has a low Sobol index so it gives only a minor contribution to the output uncertainty.
None of the strategies analyzed here satisfied our criterion of robustness with the given input distributions. This means that the probability that the number of prevalent patients in IC is larger than the IC capacity is rather high for all four exit strategies. To achieve more satisfactory performance (as far as the model and the chosen QoIs are concerned), parameter distributions corresponding to stricter interventions are needed. This would require that the effects of policy measures can “push” the bulk of the parameter distributions to more favorable values compared to the distributions used in our analysis here. In particular the insights about the safe operating space can be useful to determine how the parameters distributions might need to be modified in order to obtain more desirable outcomes, e.g. towards stricter interventions in case of Flattening the Curve or towards longer intervals between consecutive phases for Phased Opening.
The analysis presented here can be extended to a broader set of models and diseases. Since we used non-intrusive methods, the same type of analysis can be applied to a different transmission model for COVID-19 or to a computational model for a different epidemic. Furthermore, the set of methods used for analysis can be enlarged in several ways. In cases where the input parameters are not mutually independent, one can perform SA by computing the Kucherenko indices [31] instead of the Sobol indices. Also, Bayesian inference can be used to estimate the distributions of (some of) the input parameters from data, or update these distributions as more data become available [7,8]. Finally, with Bayesian model averaging [32, 33] one can address model-form uncertainties (sometimes referred to as model structure uncertainties), and therefore compare the outcome of different models for the same epidemic.
A similar model-based analysis of uncertainties has been performed by Davies et al [34], who used multiple realizations of their stochastic model in combination with variations in the basic reproduction number. Yet our analysis comprehends a larger set of uncertainties and is mathematically more rigorous as it is based on the theory and concepts of uncertainty quantification and sensitivity analysis. UQ techniques are not often applied to epidemiological models; an exception is the recent study by Edeling et al [35]. However the analysis presented there has a different scope from our study: the analyzed model and the UQ methods differ, moreover our study encompasses a broader set of strategies and addresses extensively the intrinsic model uncertainty.
It is important to realise that UQ and SA results are conditioned on the choices of parameter distributions and should therefore be interpreted with caution. As such, our results should not be interpreted as a definitive formal ranking of the analyzed exit strategies, as these strategies might show better or worse performance when considering different distributions. We aimed at picking plausible distributions, however we do not claim a homogeneous level of rigor among the choices that were made. On a related note, we assumed in this study that once the strategy is started it is not changed, implying that input parameters (or their distributions) do not change over time. Furthermore, we limited the scope of our analysis to a small number of key parameters for which we assumed specific probability distributions. However, the computational model has additional parameters (relating both to policies and to other aspects), which were kept fixed here but could be added to the set of uncertain parameters in a more comprehensive UQ analysis.
For the policy-related parameters, the chosen distributions correspond to the (assumed) effects of policy measures in the real world. The feasibility of implementing such measures is a different matter, beyond the realm of mathematical and computational modeling and therefore not considered here, but important nonetheless. As a concrete example, Fig. 5 shows that the delay in the identification of infected individuals must be short for the number of IC patients to stay below IC capacity. This result agrees well with Hellewell et al. [36], who found that when the delay between symptoms onset and isolation increases, the probability to keep the spread of the virus under control decreases. However, achieving such high levels of effectiveness may prove challenging in reality, given the increasing level of population fatigue with regard to policy adherence.
We conclude with some remarks about scaling up the analysis performed here to more extensive assessments and to more complex models (with higher computational cost). Scaling up will enable fast, frequent and comprehensive analysis of uncertainties and sensitivities in epidemiological models. Executing such analysis in a timely fashion is essential for it to be useful for policy makers. In this study we limited the analysis to a handful of uncertain parameters (less than 10), keeping all other parameters fixed. A more comprehensive study, on a larger set of parameters, will be computationally more demanding as it typically requires more model runs. Thus, access to sufficient computational resources is important to scale up.
For the analysis reported here, we had access to a supercomputer of the Poznan Supercomputing and Networking Center. A single campaign with circa 10000 runs for the SA of an exit strategy took several hours in total (including the time needed for the job submission to the supercomputer, parallel execution of the model runs on a single node with 28 cores, and retrieval of the results). If quantification of uncertainties is to be performed frequently and rapidly (e.g. in an “operational” setting with a daily or weekly cycle of producing forecasts with quantified uncertainties, or for weekly re-evaluation of a multitude of policy options), a dedicated computational infrastructure is recommendable to have uninterrupted access and to avoid long queuing times for compute jobs.
Besides access to computational resources, software suitable for efficient UQ and SA is needed. The open source VECMA toolkit used in this study is developed for use on high-performance computing platforms. Last but not least, a dedicated team with combined expertise (UQ and SA; epidemiology and computational modeling; high-performance computing and software) will be central to successful upscaling and thereby to support policy making with timely information about uncertainties and sensitivities of model results.
Conclusions
In this study we analyzed the uncertainties and sensitivities of an agent-based transmission model for the COVID-19 epidemic under four different exit strategies. Our analysis showed that the uncertainties in the model simulation results for each considered exit strategy are substantial. They were found to be mostly generated by uncertainties in the parameters directly related to the strategy itself (such as implementation and uptake of the strategy) rather than uncertainties due to other factors (such as duration of infectiousness). With the parameter distributions that we choose, the Contact Tracing strategy was the most effective. Finally, because we used non-intrusive methods, our analysis can easily be extended to other strategies as well as to other computational models and epidemics.
Data Availability
The transmission model considered in the analysis is publicly available on GitLab. The scripts used for the execution of the UQ and SA campaigns and post-processing analysis are publicly available on GitHub. All data considered in this study have been generated by the aforementioned scripts and computational model.
Supporting information
Author contributions
Conceptualization DC, FG, LEC, WE, BS
Data curation FG
Formal analysis FG
Funding acquisition DC, LEC, SJdV
Investigation FG
Methodology DC
Project administration DC
Resources DC Software FG, WE
Supervision DC
Validation FG, BS
Visualization FG
Writing - original draft FG, DC, LEC, WE, BS, SJdV
Acknowledgments
FG, WE and DC were supported by the European Union Horizon 2020 research and innovation program under grant agreement #800925 (VECMA project). The calculations were performed at the Poznan Supercomputing and Networking Center. LEC acknowledges funding from the Dutch Research Council (NWO, grant 016.Veni.178.023). LEC, SJdV and DC received support from ZonMw grant #10430022010001. The authors declare no competing interests.
S1 Appendix
Numerical implementation of the strategies
All strategies share the initial sequence of interventions which represent
the spreading of the virus among the population
the beginning of the general lockdown
the imposition of a stricter lockdown
the reopening of the schools
the application of the proposed exit strategy.
The start of the epidemic is simulated by randomly seeding 10 infections. When the amount of infected agents reaches a certain value (model parameter inc_cum_cond), a general lockdown is established. One week later a stricter lockdown is imposed (roughly mimicking the history of the Dutch lockdown in spring 2020) until schools are reopened 53 days later. Thirty days later the selected strategy is applied.
From a numerical point of view, the interventions are described by the parameters intervention_t (a vector of time points at which interventions change) and intervention_effect (a vector of the time-specific effect of interventions - i.e. a multiplier - for the average transmission rate in the population). These two parameters have to be defined for all strategies and additional strategy-dependent variables might be required (see Methods and model section in the main text).
The initial sequence of interventions shared by all strategies is implemented as follows:
intervention_t = cumsum([0, 10, 7, 53, 30, α]) where cumsum indicates that the components of the resulting vector are the cumulative sums and the second component is dynamically adapted such that the lockdown starts when there are inc_cum_cond infected individuals;
intervention_effect = [1, 0.3, 0.15, 0.25, β].
We use the Greek letters α and β to indicate either single values or sequences of values that define the strategies. More details about their specific definition in the cases considered follow below. α is necessary only in case of Intermittent Lockdown or Phased Opening and can be omitted in Flattening the Curve or in Contact Tracing.
Flattening the Curve
For this strategy it is not necessary to extend further the time series of interventions, i.e. intervention_t = cumsum([0, 10, 7, 53, 30]), and β is modelled by the uncertain parameter int_effect such that intervention_effect = [1, 0.3, 0.15, 0.25, int_effect]. We assume that the population uptake is constant in time, hence uptake is kept unchanged for the whole simulation, i.e. intervention_uptake = rep(uptake, len(intervention_t)) where uptake is repeated as many times as the length of intervention_t
Contact Tracing
Also in this strategy we do not require α and the time series of interventions is given by intervention_t = cumsum([0, 10, 7, 53, 30]). Since the restrictions imposed by CT are determined by the variables trace_rate_I, trace_prob_E and trace_contact_reduction and not by intervention_effect, we set β = 1, i.e. intervention_effect = [1, 0.3, 0.15, 0.25, 1]. We assume that the CT measures are not active before the actual start of the strategy, hence we have trace_prob_E = [0, 0, 0, 0, trace_E], trace_rate_I = [0, 0, 0, 0, trace_I] and trace_contact_reduction = [0, 0, 0, 0, trace_cr] where trace_E, trace_I and trace_cr are the uncertain parameters discussed in the main text.
Intermittent Lockdown
The interchange between lockdowns and lifts is modelled by intervention_t and intervention_effect. In particular, α consists in a phase of lift after the opening of the schools - necessary to let the virus spread again among the population - and the repetition of the length of the lockdown and subsequent lift. The duration of the very first lift phase after the general lockdown is largely determined by the initial state of the epidemic at the start of the strategy, thus it is not useful to consider it as uncertain and we set it to 25 days for our analyses. The final time series of the interventions reads intervention_t = cumsum([0, 10, 7, 53, 30, 25, rep((lockdown_length, lift_length),10)]) where (lockdown_length, lift_length) is being repeated 10 times.
Similarly β consists in the alternation of the effect on the average contact rate due to the lockdown or the opening, i.e. intervention_effect = [1, 0.3, 0.15, 0.25, rep((1, lockdown_effect),11)] where (1, lockdown_effect) is repeated 11 times to match the length of intervention_t. As for Flattening the Curve, uptake is assumed not to vary during the simulation, i.e. intervention_uptake = rep(uptake, len(intervention_t)).
Phased Opening
The phased lift plan is generated by the function gen_phased_lift and does not require extra ad-hoc numerical coding. This function includes several input parameters, of which only some have been considered here as uncertain, i.e. pl_intervention_effect_hi, intervention_lift_interval and uptake. We refer the reader to the documentation of the virsim model provided in its GitLab repository [6] (script gen_intervention.r in the R folder) for a detailed description of the gen_phased_lift function, its inputs and its outputs.
S2 Appendix
Theory of Sobol indices
Here we will briefly outline the theory and the algorithm by which we compute the Sobol indices. For a more detailed description we refer to [19]. Let Y = ℳ (X1,, Xp) be the output of a computational model, with p independent inputs such that the joint pdf p(X) is given by Then, the output variance, conditional on a fixed value (where j ∈ {1, …, p}), is written as Note that Ω∼j indicates integration over the support of all inputs pdfs, except that of p(Xj). Similarly, Ωj denotes integration over only the support of p(Xj), and we will use Ω to denote integration over the entire stochastic input domain. Also, similar to the main text, we write 𝔼i and 𝕍ari for expectation and variance over Xi only, 𝔼∼i and 𝕍ar∼i for expectation and variance over all X1, X2, …, Xp except Xi, and 𝔼 and 𝕍ar for expectation and variance over all X1, X2, …, Xp.
Expression (1) allows one to gauge the effect that keeping one input fixed has on the output variance. Sobol indices, defined later on, are global sensitivity measures however, so we wish to eliminate the dependence on the specific, local value . We therefore integrate the conditional variance (1) over all values, i.e. Again, this is the expected value of the variance of Y, while keeping one input (Xj) fixed. Hence, if we subtract this from the total variance of Y, we get the contribution to the variance due to Xj alone: where 𝕍ar[Y] := 𝔼[Y 2] 𝔼 2 [Y]. We can write the above difference between the variance and the expectation of the conditional variance, as the variance of the conditional expectation: which equals (2). The first-order Sobol indices Sj are defined as the ratio of the variance due to Xj alone, i.e. (3), over the total variance. This gives Note that this is the expression found in the main text. The appearance of the term 𝕍arj [𝔼∼j [Y | Xj]] might (suggest an implementation involving a double loop, where the inner loop computes , and the outer loop computes the outer integral. If we denote M as the number of MC samples, the involved cost would be M 2. However, Saltelli [19] developed an algorithm, through which the cost of computing the first-order indices is reduced to M (p + 1). Briefly, several input matrices are constructed, of which the first two are filled with independent draws from p (X): We can think of A as the ‘sample’ matrix, and B as the ‘resample’ matrix. Next, p matrices are constructed, where the j-th column of B is replaced by the j-th column vector of A: We can now approximate the first-order Sobol indices as follows: where This single-loop MC approximation corresponds to the integral , see [19] for details.
We can estimate the unconditional mean and variance in (6) by evaluating the model ℳ on the M samples of A. To estimate the Uj, we can see from (7) that this requires evaluating the model on the rows of Bj, j = 1, …, p, where all inputs except Xj are resampled. This brings the total cost up to M (1 + p). If one wishes to also compute the total-order indices, the cost will be increased to M (p + 2), see [19] for details.
S1 Fig. Distributions of the uncertain input parameters.
S2 Fig. Total Sobol indices of the first QoI.