Optimal algorithms for controlling infectious diseases in real time using noisy infection data

Optimal epidemic control

Our study focuses on epidemic control at the early stages of an epidemic with no or limited immunity in the population and without any available pharmaceutical remedies. Consequently, the main control measures are NPIs, such as social distancing, mask-wearing, stay-at-home orders, business closures and travel restrictions. These measures limit disease transmission with the aim of reducing the likely numbers of severe infections below healthcare capacities and minimising expected morbidity and mortality. However, these benefits must be balanced against the costs induced by those NPIs, which may include economic downturns and loss of livelihood.

An ideal control policy would keep the incidence of new infections at a manageable (target) level, optimally balancing the costs of treating infections and implementing interventions. However, this is non-trivial both because disease dynamics can change in real time and our ability to track those changes is strongly limited by the quality of available surveillance data. To achieve this, we propose a model predictive control (MPC) framework for optimising epidemic interventions based on real-time incidence data. MPC utilises a mathematical model to project the dynamical behaviour of the controlled system [29].

We outline our MPC approach in panel (a) of Fig. 1. Broadly, this algorithm compares historical incidence data to desired target levels and decides appropriate control actions that minimise both expected future costs and drive us closer to our target. This target is essentially the level of infections policy-makers may believe sustainable. For example, this may be set so that the proportion of severe infections leading to hospitalisation never exceeds healthcare capacities or to ensure a manageable endemic level of infections. This joint target-cost optimisation process is iteratively done in real-time via the feedback loop in Fig. 1 and makes use of short-term projections.

• Download figure
• Open in new tab

Figure 1:

Panel (a): Schematic diagram of model predictive control for optimising epidemic interventions. The top chart introduces the elements of the feedback loop where the actions of the agent are chosen according to the incidence of new infections (or a proxy such as new cases) which is the monitored output state. The highlighted panel explains the model predictive method for selecting the optimal NPI from the action space which is based on using short-horizon projections and the expected reward, E(ρ), under those projections for various strategies. The strategy with the highest reward (minimum cost) is implemented. The reward or cost here usually depends on how far the epidemic state is from our desired objectives, which may aim at jointly reducing both severe epidemic outcomes and intervention intensity. Panels (b) and (c) show event-triggered and time-triggered alternative strategies, with NPIs implemented or relaxed based on incidence thresholds (the event trigger) or according to a pre-defined schedule (the time trigger). Panel (d) illustrates how realistic surveillance imperfections such as reporting delay and under-reporting distort the true incidence of infections into the incidence of cases, which we practically must use to inform decision-making. Blue curves represent the reported cases, while red curves indicate true infection incidence. Reporting delay manifests (approximately) as a time-lag with respect to the true incidence curve, whereas under-reporting results in a stochastic downscale of the incidence curve along the vertical axis.

Our algorithm is analogous to a Markov Decision Process [30] and involves an agent (policy-maker) that decides which NPI to implement from a finite action-space of possible NPIs, based on the projected discounted reward. This reward is compounded over a fixed time horizon for each possible action or action-sequence, as explained in the highlighted incidence curve in panel (a) of Fig. 1. In our simulations, we consider three possible actions which we will refer to as full lockdown, limited social distancing and no restrictions. These interventions are modelled as multiplicative reductions in transmission. The reward function is derived from the costs associated with infections and interventions and is evaluated using projections informed by daily new cases data (C_t) and the expected changes in transmission due to our possible intervention choices.

We use cases as we rarely know the actual number of infections I_t due to surveillance imperfections such as under-reporting and delays. The overall cost consists of the cost attributed to the interventions, a penalty term that is proportional to the error between the incidence and our target and an additional penalty term applied for overshoots larger than 50% of the incidence target. The intervention cost is accumulated daily, with large, moderate and no cost for each day under full lockdown, limited social distancing and no restrictions, respectively. At each decision point, the action with the largest projected discounted reward (i.e., the smallest cost) is chosen. We only consider scenarios in which these decision points occur at weekly or larger spacings to model practical review periods.

We focus on three key metrics to evaluate the performance of control algorithms informed by epidemic data that is subject to different sources of surveillance noise. We look at the peak (maximum) of the incidence curve, the bounding envelope (difference between the maximum and minimum) of incidence values in periods when the epidemic is stable (akin to the endemic state) and the cost spent on interventions. According to these metrics a high-performing control algorithm prevents large epidemic peaks, and forces infections into a manageable steady state with small fluctuations without overly applying NPIs. The last is achieved by optimising timing, switching and duration within our available NPI set. The mathematical formulation of our MPC approach is presented in the Methods.

Alternative control strategies

We compare the performance of the above proposed optimal control algorithm with two simpler control strategies: an event-triggered feedback control and a cyclic time-triggered control strategy. These strategies represent two fundamental approaches for controlling epidemics in real time, which have either been applied or proposed earlier. Note that for all strategies that we consider, we allow tuning so that the strategies can stabilise the observed incidence.

The event-triggered control applies or relaxes lockdowns whenever reported incidence crosses a predefined threshold, which is often heuristically set in practice. This crossing constitutes an event. We illustrate this approach in panel (b) of Fig. 1 for an example incidence time series. Event-triggered control approaches have previously been used to enact interventions, e.g. for influenza [31] and were considered for triggering NPIs to suppress COVID-19 in the UK [3]. Although this strategy applies limited feedback based on the most recent incidence, it is unable to leverage the information in the full past time series or assess the likely future outcomes of its decisions.

In contrast, the cyclic time-triggered control strategy (see panel (c) of Fig. 1) implements a predefined sequence of actions, which is not based on any direct feedback from the epidemic dynamics. The periods with full lockdown or no restrictions can be arbitrarily long, e.g. a 20/10 cyclic control strategy repeats 20 days of full lockdown followed by 10 days of no restrictions. This strategy was proposed as an effective means of COVID-19 control when surveillance data are poor quality and hence unreliable for informing decisions [24].

Observation noise and uncertainty

Infection data are rarely observed directly. As a result, a proxy such as the incidence of reported cases is commonly used for informing decisions in real time. Cases can be modelled as noisy versions of infections [32, 33] that are subject to under-reporting and delays in reporting. We examine how the performance of the above control strategies varies with the levels of these noise sources.

Reporting delay describes the time between an infection and when it is reported as a case. This subsumes lags between infection and symptom onset, onset and presentation at a health facility and processing time to confirmation (e.g., via testing). As explained in panel (d) of Fig. 1, the reporting delay is perceived as a probabilistic shift in time of the reported cases when compared to the actual infection data.

Under-reporting or under-ascertainment occurs when not all infected individuals are tested and reported as cases. This may result, for example, when infections are asymptomatic or only cause mild symptoms so that infected individuals do not seek medical care or testing. Under-reporting can be modelled as a (stochastic) downscale of the infections i.e., a random fraction of new infections appear as new cases.

In the following sections we explore optimal MPC, event and time triggered control strategies and their performance as under-reporting and delays increase. Full details (including mathematical descriptions) of the algorithms we employ are available in the Methods.

Optimal MPC performance

We demonstrate the performance of our MPC algorithm first on perfectly observed incidence data and then explore the influence of practical surveillance limitations. Figure 2 presents four scenarios using MPC to mitigate an infectious disease with generation time (mean: 6.5 days [34, 35]) and basic reproduction number (R₀ = 3.5 [36]) chosen to match those previously estimated for COVID-19. Row (a) shows the ideal case where the agent has access to the true incidence of infections. Here, the MPC strategy is able to keep incidence near the target incidence level (5000 new infections/day). Note that we cannot precisely track the target even in this ideal data setting due to intrinsic stochasticity in the epidemic, a finite action space of possible interventions and a policy review period that is at least a week. The resulting fluctuations about the target are a measure of the fundamental control performance under these settings.

• Download figure
• Open in new tab

Figure 2:

Simulation results with optimal epidemic control. The left panels show reported cases from ensembles of 100 simulations using generation times and reproduction numbers estimated for COVID-19. The faded curves show different individual realisations of the epidemic with the three black curves marking the 5% and 95% percentiles of the ensemble and the mean reported cases. The horizontal dashed line shows the incidence target. The highlighted thick curve of reported cases is coloured based on the NPI implemented on a given day. The coloured thin curve indicates the true incidence corresponding to that highlighted simulation. The right column shows similar diagrams for the effective reproduction number. In that column the faded grey curves and the highlighted curve represents the estimated effective reproduction number from the true infection incidence whereas the thin curve indicates the estimated value from the reported cases in the highlighted realisation of the epidemic. The inset pie charts in the right column indicate the ratio of days spent under a given NPI across the full simulation ensemble. Row (a) represents the baseline case without delay or under-reporting with a policy review period of t_rev = 7 days. Compared to the baseline case, panel (b) shows simulations with reporting delay (mean 7 days, shape parameter α = 5.0), panel (c) simulations with under-reporting (with mean reporting rate 0.3, dispersion a = 8.0). Panel (d) has no observation noise but provides simulations under an increased policy review period of t_rev = 14 days.

In rows (b) and (c) we demonstrate how reporting delay and under-reporting separately affects MPC (cf. panel (d) in Fig. 1). As indicated by the panels in the right column, surveillance imperfections of infection incidence also show up in the estimated reproduction number leading to discrepancies between the observation and the true state of the epidemics. Row (b) shows the effect of reporting delay with a mean of 7 days. In this case, the MPC strategy is still able to keep incidence at a manageable level, however, the delay in observing cases results in late intervention that spurs a higher peak in incidence and larger fluctuations once the epidemic is under control. The thick highlighted curve shows the cases that inform the agent or decision-maker while the thin highlighted curve are the true (unknown) infections.

Row (c) shows the effect of under-reporting with a mean reporting rate of 0.3, i.e., only 30% of infections are reported as cases. Our MPC algorithm is able to achieve the target level but the true incidence fluctuates at a higher level due to the under-reporting noise process which, on average scales down infections to reported cases according to the mean reporting rate. The stochasticity of the under-reporting also occasionally misleads the MPC strategy to believe that the epidemic is under control and the incidence is below the target level or the effective reproduction number is smaller than its true value, causing higher epidemic peaks than in the ideal surveillance scenario. Interestingly, the average proportion of NPIs across the simulation time horizon is similar across scenarios (a-c), indicating that the deviations in the performance of the optimal MPC strategy are due to mistimed interventions resulting from the noise in surveillance, causing either overly conservative or relaxed policies.

Row (d) in Fig. 2 shows the effect of increasing the policy review period from 7 to 14 days with no observation noise. In this case, the MPC strategy is still able to keep incidence at a manageable level, however, the fluctuations in daily incidence, and consequently the peak and the bounding envelope of the later stabilised epidemic are larger as the agent has a reduced ability to intervene and update control actions. The overall effect of increasing the policy review period is similar to having a reporting delay as ultimately, both lead to delayed responses.

The limits of control due to delayed reporting

We assess how the delay in reporting limits the performance of each control strategy. We consider scenarios with different but stationary reporting delay distributions. The delay for a single infection follows a Gamma distribution τ ∼ Gamma(α_τ, β_τ), with shape and scale factors α_τ and β_τ, respectively. The mean reporting delay is then α_τ β_τ while the variance is. This is a common model of reporting delays and has been used to describe surveillance of COVID-19 and Ebola virus disease among others [34, 37, 38, 39]. To control the mean delay τ_mean and dispersion α directly we re-parametrise the distribution by the choice α_τ = α and β_τ = τ_mean/α. This means that for a given mean reporting delay, the variance is inversely proportional to the dispersion parameter α, i.e., larger values of α correspond to more deterministic reporting delays.

We consider 6 mean reporting delays ranging from 3.5 to 21 days, each with 4 different variances. This allows us to characterise how the mean and the dispersion of the reporting delay limit optimised and heuristic control strategies. For each scenario, we run an ensemble of 1000 simulations, each with a different random seed. Fig. 3 plots key results for simulations under COVID-19 disease parameters. Each column of panels depicts the performance of a different control strategy: MPC, event-triggered, and time-triggered cyclic control. The first row plots the peak incidence, the second row shows the size of the steady-state solution envelope (the range of oscillations after the epidemic is under control) and the third row charts the average costs of NPIs. When calculating the steady-state envelope, we look at the maximum and minimum of daily cases after the incidence of cases falls below the target and the effective reproduction number is below 1. It is not possible to determine a steady-state envelope for every epidemic realisation of every scenario as under some parameter combinations the control algorithm may fail to stabilise the outbreak. If this happens we use the peak value as the size of the solution envelope. We find that the distribution of peaks and steady-state envelopes are multi-modal for certain parameter combinations. This results from applying a controller that has a discrete action space and policy review time longer than the algorithm time step (which is daily).

• Download figure
• Open in new tab

Figure 3:

The impact of increasing reporting delay on optimal control. The left panels show results for our MPC algorithm, while the middle and right panels respectively present equivalent outputs under event-triggered (lockdown and relaxation thresholds of C_LD = 2500 new cases/day, and C_relax = 1500 new cases/day) and time-triggered (cycle of 45 days lockdown, 9 days of no restrictions, starting on day 38) strategies. Row (a) shows scatterplots for peak incidence and row (b) illustrates the steady-state envelope size relative to the target incidence across different time delay distributions. The horizontal axis represents different mean reporting delays with colours depicting different dispersion levels. Larger values of the parameter α indicate more deterministic delays. Row (c) shows the mean intervention costs for each ensemble of 1000 epidemics, simulated under estimated parameters from COVID-19. An equivalent analysis for simulated Ebola virus disease epidemics is presented in Fig. 7 of the Supplement.

For the comparison in Fig. 3, we tuned the parameters of the different control algorithms such that they involve similar intervention costs. Additionally, we selected a rolling policy of 45 days of full lockdown followed by 9 days of no intervention (i.e. a 45/9 cyclic policy), which is enacted 39 days after the simulated epidemic started. This is to ensure that the long term average effective reproduction number is approximately 1, so we observe the epidemic reaching a near steady-state. Arbitrary choices of policy lengths may fail to stabilise the epidemic, resulting in continued (but slower) growth or rapid suppression that is costly.

Figure 3 indicates that reporting delays have a detrimental effect to feedback-control strategies i.e., the MPC and event-triggered control approaches, as mistimed action leads to drastically higher incidence peaks and wider steady-state envelopes. In some realisations, the desired target is exceeded by a factor as large as 100. Interestingly and perhaps counter to intuition, the detrimental effect of reporting delay is less if the variance of the reporting delay is high. Low variance means that the reporting delay is almost deterministic (see results with higher values of the dispersion parameter α), while high variance implies that we have a larger portion of cases where the delay is small (and also more cases subject to very large delays, for the same overall mean lag). As a result, high variance delay distributions allow some access to more recent infection information, which can aid decision making. Row (c) of Fig. 3 shows the average costs spent on interventions for each strategy under different delay settings. As expected, we find the optimal control strategy is notably more economical than either the event-triggered or the time-triggered method.

Time-triggered cyclic control is insensitive to the reporting delay as it follows a pre-defined sequence of actions irrespective of the observed incidence trajectory. As a result, while the costs of the interventions are higher than those of MPC, it is possible to devise predefined intervention strategies that outperform optimal feedback-control under large case reporting delays with means of 10.5-14.0 days or more, depending on the delay dispersion. This exposes the limits that surveillance quality can impose on our ability to reactively control an epidemic. Importantly, as long as the dispersion of our reporting delays is sufficiently large (α ≈ 1) then optimal MPC strategies offer a cost-effective intervention approach for any mean delay.

In the Supplement, we also present the results of the same analysis for Ebola virus (See Fig. 7). We find that the performance deterioration caused by reporting delays follows similar trends to what we observed for COVID-19 in Fig. 3. However, a key difference is that Ebola virus has a longer generation time, which results in slower dynamics, that reduce some of the negative impact of reporting delays. For comparison, with the MPC approach and a near deterministic (α = 200) delay distribution with a mean of 21 days, we obtain peak infection incidence 100-1000 times larger than target for COVID-19. In contrast, the same delay only leads to peaks of about 4-5 times the target for Ebola virus disease.

The limits of control due to under-reporting

Having isolated the influence of reporting delays, we now examine the impact of under-reporting. This is another common and important surveillance imperfection in which not all newly infected individuals are reported as cases leading to under-ascertainment of the true numbers of infections and hence the size of the epidemic. We model the number of daily reported cases with a Beta-binomial distribution C_t ∼ BetaBin(I_t, α_ν, β_ν), with shape parameters of α_ν and β_ν. The expected number of reported cases is then I_tα_ν /(α_ν +β_ν) while the variance is I_tα_ν β_ν (α_ν +β_ν +I_t)/((α_ν +β_ν)²(α_ν + β_ν + 1)). We refer to the ratio of the expected reported cases and the true infection incidence as the mean reporting ratio ν := E[C_t/I_t] which is constant in time. To directly control the mean reporting ratio ν_mean and the dispersion, we choose α_nu = a and β_nu = a(1 − ν_mean)/ν_mean. This means that larger values of the dispersion parameter a result in a smaller dispersion in reporting rate.

We consider 6 different mean reporting ratios from 0.1 to 0.85 each with 4 different variances. A smaller reporting fraction means larger under-reporting. Fig. 4 plots the performance of the MPC, event-triggered and cyclic strategies in successive columns. For these diagrams, the target incidence is scaled up according to the mean reporting ratio to have comparable results.

• Download figure
• Open in new tab

Figure 4:

The impact of increasing under-reporting on optimal control. The left panels show results for our MPC algorithm, while the middle and right panels respectively present equivalent outputs under event-triggered (lockdown and relaxation thresholds of C_LD = 2500 new cases/day, and C_relax = 1500 new cases/day) and time-triggered (cycle of 45 days lockdown, 9 days of no restrictions, starting on day 38) strategies. Row (a) shows scatterplots for peak incidence, row (b) illustrates the steady-state envelope size relative to the target incidence for different reporting rate distributions. The horizontal axis represents different mean reporting ratios with colours depicting different dispersion levels. Larger values of the parameter a belong to more deterministic case reporting distributions (i.e., constant reporting). Row (c) shows the mean intervention costs for each ensemble of 1000 epidemics, simulated under estimated parameters from COVID-19. An equivalent analysis for simulated Ebola virus disease epidemics is presented in Fig. 8 of the Supplement.

We find that as the variance in under-reporting increases the performance of optimal feedback strategies deteriorates. This follows because the under-reported case curve is effectively a stochastic downscaling of the true infection incidence curve. As a result, larger stochasticity in fluctuations for a given mean reporting rate more substantially distorts the observed incidence and misleads control. Low variance means that the under-reporting is more deterministic, so that reported case incidence better resembles the shape of the true infection incidence curve.

Both our MPC and the event-triggered strategy show similar patterns in how the peak and envelope of the controlled epidemics vary with the under-reporting statistics. However, the MPC algorithm achieves better performance with smaller intervention costs. This improvement derives from the the MPC approach leveraging all the available historical information in the incidence curves.

Time-triggered cyclic control is not affected by under-reporting because it is agnostic to the incidence data and how it is reported. The apparent variation in the performance of the cyclic control strategy is due to the scaling of the target according to the under-reporting rate. While under moderate noise and uncertainty, our MPC approach is more effective and cost efficient than either of the heuristic methods, cyclic control can outperform MPC if the reporting rate is very low and the variance is high.

Comparing the effect of under-reporting on control of COVID-19 and Ebola virus (see Fig. 8 in the Supplement) we find similar trends for both pathogens. Even though the drop in performance that we observe is less pronounced than that due to reporting delay, the slower dynamics of Ebola virus still makes it more controllable than COVID-19 when uncertainty in case under-reporting is present. This is evidenced by the solution peaks and steady-state envelopes being distributed in a smaller interval across the simulation ensembles for Ebola virus disease than for COVID-19.

Integrating noise and intervention frequency

Having explored the performance limits induced by reporting delays and under-reporting in isolation, we now consider their combined impact. We analyse ensembles of simulations under realistic noise distributions for COVID-19 and Ebola virus disease. The uncertainty in case reporting data can vary markedly depending on the context and national or regional differences in how surveillance is conducted. In [1], reporting delays of 9-12 days were estimated for COVID-19 in Italy, whereas [40] inferred case-reporting rates between 7-38 % across France. Based, on these, we consider a mean reporting delay of 10.5 days with a dispersion of (α = 5.0) and a mean reporting rate of 0.3 with a dispersion of a = 8.0.

Rows (a) and (b) in Fig. 5 show results for simulated COVID-19 epidemics under realistic noise distributions and subject to policy review periods of 7 and 14 days, respectively. We find that the MPC strategy stabilises the epidemic around the target, but fluctuations are considerably larger than in the baseline case (see Fig. 2a) due to the delay and under-reporting. This results in an overshoot of about 5-times the target incidence under a 7-day policy review period and around 10-times for a 14-day policy review period. Comparing these results with the simulations in panels (c) and (d) of Fig. 5 for epidemics simulated under Ebola virus parameters, we observe that, the MPC strategy is more effective in controlling these epidemics and the effect of noise is less detrimental to controllability. For the Ebola virus epidemics (panel (c)) we find that running the MPC algorithm with a 7-day policy review period causes a 3-times overshoot of the target incidence of cases, which remarkably, only slightly deteriorates to about 4-times when a 14 policy review period is used (see panel (d)). This occurs because the longer generation time of Ebola virus results in a slower epidemic that allows the MPC algorithm more time to effectively adapt to changes in the epidemic dynamics. Consequently, we must act more swiftly when responding to diseases with shorter generation times as their faster growth can quickly destabilise data-informed policies. Note that in scenarios where the same epidemic parameters were used ((a) and (b) – COVID-19, (c) and (d) – Ebola virus disease) we observe that our MPC algorithm enacted and sustained NPIs in similar time ratios (with COVID requiring more restrictions than Ebola). This confirms that differences in performance are due to timing instead of overly conservative or relaxed strategies.

• Download figure
• Open in new tab

Figure 5:

Optimal control and policy review for realistic simulated epidemics. Rows (a) and (b) represent epidemics simulated under a COVID-like generation time and basic reproduction number with 7 and 14 days policy review periods, respectively. Rows (c) and (d) represent epidemics simulated under Ebola-like parameters with 7 and 14 days policy review periods, respectively. Because Ebola virus disease has a longer generation time, we discard a burn in period of 14 weeks to allow the epidemic to grow towards the initial target. The left column show reported cases from ensembles of 100 simulations with mean delay 10.5 days, delay dispersion α = 5.0, mean reporting rate 0.25, a = 8.0. These settings reflect estimates of realistic surveillance noise from the literature. The faded curves show different individual realisations of the epidemic with the 3 black curves marking the 5% and 95% percentiles of the ensemble and the mean reported cases. The horizontal dashed line shows the incidence target. The highlighted thick curve of reported cases is coloured based on the NPI implemented on a given day. The coloured thin curve indicates the true incidence corresponding to that highlighted simulation. The right column shows similar diagrams for the effective reproduction number. The faded grey curves and the highlighted curve represents the effective reproduction number estimated from true incidence whereas the thin curve indicates the estimated value from the reported cases for the highlighted realisation of the epidemic. The inset pie charts in the right column indicate the ratio of days spent under a given NPI across the full simulation ensemble.

Epidemic governing equations

We model the spread of the disease in a population as a generalisation of the standard renewal branching process [28]. This model is used both to make projections that inform optimal control and to simulate ‘ground truth’ epidemic trajectories. The renewal branching process is a stochastic model describing how the incidence of new infections on day t, I_t, depends on past infections at times s ≤ t and the characteristics of the disease. This is captured by the Poisson distribution where R_t is the effective reproduction number on day t, with the set of weights, w_t−s for all s, obtained from the generation time distribution [34] of the disease. We assume that the generation time distribution is known or estimated from other paired transmission data. The weight w_t−s is the probability that a secondary infection occurs t − s days after its primary infection. As is standard practice, we model the stochasticity of the generation time with a Gamma distribution

The shape and scale factors α_gen and β_gen parametrise the probability density function p_gen(t). The weights w_t−s used in Eq. (1) are then calculated as. We consider two generation time distributions, with respective parameters provided in Table 1, which are commonly used to describe epidemics of Ebola virus disease [51] and COVID-19 [35].

The above formulation generalise s the standard renewal model, which describes incidence I_t ∼ Pois with the sum referred to as total infectiousness of all infectious individuals [52]. However, this classical formula implies that any intervention applied to curb the spread of the disease results in an immediate change in the reproduction number. We apply the generalised formula of Eq. (1) to model scenarios where the infectiousness of a population and the impact of interventions depend on both the history of infections and reproduction numbers. The latter dependence has a smoothing effect that allows for the finite time effects of realistic interventions on transmissibility. A similar generalisation was introduced in [33].

The effective reproduction number is derived from the basic reproduction number R₀, which describes how many people an infected individual is expected to infect in a fully susceptible population. When an NPI is introduced the basic reproduction number R₀ is multiplicatively changed by a factor c_t, yielding

We consider an action space with three possible NPIs: no intervention (c_t = 1), limited social distancing (c_t = 0.5) and full lockdown (c_t = 0.2). While this is a simplified classification of NPI types, our control framework allows for more possible actions to be easily modelled. Although in reality the factor c_t is unknown and needs to be estimated, we assume throughout this study, that the effect of any NPI on the reproduction number is known and without uncertainty. Our framework does allow the inclusion of uncertainty on these effects and we present analyses under stochastically varying c_t in the Supplement. However, our goal is not to describe precisely how interventions attenuate transmissibility, but instead to derive insights into how noise influences the optimal timing of interventions generally. Thus, our choices of c_t are meant to be sensible but not exact and we generally do not include the uncertainty on c_t to isolate the influence of the surveillance noise. Note that if estimates of c_t are available or derived from auxiliary data, these can be seamlessly integrated within our framework for more precise results.

Optimal model predictive control (MPC) of epidemics

The model predictive control algorithm we propose aims to curb disease spread, jointly minimising the risks and costs arising from infections and applied interventions. We outline our control framework in panel (a) of Fig. 1, which consists of the following elements: a plant (the controlled system) with observable states, state-transition probabilities, an agent with an action space defining possible control actions, and a reward function.

In our model, the plant is the population where the disease is spreading, while the output state monitored is the incidence (number of daily new infections) I_t. The state transition probabilities, i.e. the probability of transitioning to any I_t+1 from any given I_t are not explicitly defined, but are implicitly determined by the Poisson distribution of the renewal model (see Eq. (1)). The control framework we use here largely overlaps with Markov decision processes (MDPs) [53]. However, the renewal model utilises both the immediate and past incidence. This is not exactly Markov but may be reconfigured into an MDP if higher dimensional state spaces are used [54, 55].

The agent in the context of an epidemic is the public health policy-maker i.e., the individual or group responsible for proposing or removing NPIs, while the action-space comprises the possible NPI choices. We consider 3 levels of interventions that we class as no intervention, social distancing and full lockdown. This broadly models stepped interventions which were common across the COVID-19 pandemic. These include the three tier system that England used to enforce localised NPIs in 2020, the 4-level alert system applied by New Zealand and related policies taken by Italy, France, Canada and others [56, 57]. Our framework computes decisions based on the projected reward over a fixed time-horizon which incorporates the costs of possible actions in our decision space and their risks in terms of expected infections.

In the reward function, we account for costs arising from the economic impact of NPIs and the risks associated with high incidence. We consider a target incidence level that defines some manageable infection level and define the absolute error I_err = |I_t − I_target |. This target may relate to healthcare capacities e.g., setting a level of incidence such that the expected hospitalisations resulting from that incidence do not overwhelm healthcare resources. Although, studies rarely consider a target incidence level, our aim is to understand and characterise the intervention tradeoffs (e.g., timing choices) that can jointly limit expected infections and the costs of those interventions. Setting I_target = 0 is analogous to an elimination target, which models the broad aims of pandemic policies employed by New Zealand [58] and China [59], for example. An I_target > 0 recognises that elimination is difficult, particularly in the face of infection reintroductions and so refocuses on stabilising healthcare burdens to sustainable levels that balance the supply and demand of health resources. Additionally, as we want to minimise the risk of large infection peaks and overshoots, our reward function also includes a penalty term ϕ_over that activates when I_t > 1.5I_target but is zero otherwise.

While regulating disease spread within the limits of healthcare capacities is of paramount importance, interventions that restrict mobility or close businesses and trade generate substantial economic and other costs. We model this with a term ϕ_t attached to every element of the action-space. There is no cost under no restrictions and the cost of full lockdown is assumed to be 15-times larger than that of limited social distancing. While some studies into COVID-19 NPIs suggest stringent interventions are 5-6 times more costly than more limited measures [7], our factor was chosen to more markedly distinguish between our two NPI tiers so that general qualitative insights could be better derived. Including all the above components, the reward function on day t is calculated as the negative quantity

The agent’s task is to choose the action which maximises the expected reward. However, there are practical limitations to decision-making. While we use daily incidence data to inform our epidemic model, we allow policy review to only occur every 7 or 14 days, i.e., the agent can only change control actions with this frequency. The time between policy updates is t_rev and reflects practical intervention constraints, e.g., both in terms of logistics and ensuring compliance, policymakers may not want to switch NPIs any faster than weekly. We also impose a practical constraint on reward optimisation by considering only finite time horizons for assessing the costs of any action. We denote this projection horizon t_proj. This models the fact that only short-term forecasts are known to be reliable for epidemic decision-making [60].

The agent calculates the expected reward for each action by simulating the epidemic with all possible control states until t_proj and taking the total temporally discounted reward where γ < 1 is the temporal discount factor. A higher γ means that the agent is more concerned about long-term rewards, whereas a smaller γ means that shorter term benefits are emphasised. Since the reward function is stochastic, the expected reward ρ is probabilistic. We therefore compute the expected total reward E(ρ) over an ensemble of simulations. This also allows us to factor in the intrinsic variability of the epidemic generating process (e.g., from the random times between infections).

If the projection horizon is longer than the policy review period, i.e., t_pred > t_rew, then, we can also propose sequences of actions over the projection horizon to be taken by the agent. We may then compute the expected reward for each action sequence then implement the first action of the sequence with the best projected reward. However, for the scenarios we consider, this approach increases computational complexity but does not improve performance. Consequently, for these longer projections, we maintain our original approach of only evaluating single possible actions and their consequences across the horizon. We collect all the key epidemic and control algorithm parameters in Table 1.

Surveillance noise and uncertainty in incidence data

Ideally, the agent would make decisions about possible NPIs based on the infection incidence in the population. Unfortunately, infection data are rarely available and a proxy such as the incidence of confirmed cases or deaths is commonly used. We focus here on the daily incidence of cases C_t but note that other proxies have analogous descriptions [32, 33]. These proxies are commonly subject to practical surveillance imperfections, which we define via a stochastic reporting delay τ and a stochastic reporting rate ν. In our model of the surveillance imperfections, the true infection incidence data I_t is first distorted by delay, then we consider under-ascertainment of the delayed cases (see Fig. 6).

Reporting delay describes the lag between an infection and its proxy. For cases this includes latencies such as the time taken between infection and presenting symptoms or confirmation via testing. In our framework, we model the reporting delay for a single case using a Gamma distribution with shape and scale factors α_τ and β_τ, respectively. The mean reporting delay is then α_τ β_τ while the variance is. To control the mean delay τ_mean and dispersion α directly we re-parametrise the distribution by the choice α_τ = α and β_τ = τ_mean/α. The cases reported with delay on day t then result from a weighted sum of past incidence at day s and the probability that it takes t − s time units before those infections are reported or confirmed as cases. This follows as where the weight factors are derived from the Gamma distribution of the reporting delay,.

The reporting rate models incomplete or under-reporting, which captures the fact that proxies commonly represent only a fraction of infections. For example, asymptomatic and less severe infections are unlikely to appear as cases. This means that only a fraction ν of the delayed infection incidence is reported with ν_t ∈ [0, 1]. In our model, we assume that the number of reported cases follows a Beta-binomial distribution

Consequently, the expected number of reported cases is while the variance is. We refer to the ratio of the expected reported cases and the true infection incidence as the mean reporting ratio which is constant in time. To directly control the mean reporting ratio ν_mean and the dispersion a we choose α_nu = a and β_nu = a(1 − ν_mean)/ν_mean.

In some cases, we investigate the isolated effect of reporting delay or under-reporting, which means that for these simulations the other source of surveillance imperfections is turned off. If there is no reporting delay, and similarly, if there is no under-reporting. These stochastic delay [61] and under-reporting [62] models have been widely used to describe surveillance noise in the literature, as well as serve as the starting point for deconvolution and nowcasting methods that attempt to correct for these noise sources [33, 63, 64, 65].

Estimation of the reproduction number

When projecting likely infections (or proxies) over a horizon in our algorithm, we assume knowledge of the effect of NPIs on the reproduction number. This is captured by the coefficients c_t. However, we do not assume knowledge of the true basic reproduction number and so must estimate this quantity from past data.

We start by inferring the time-varying effective reproduction number by applying the formula [52] where Λ is the total infectiousness and is calculated as with weights w_s derived from the generation time distribution. Then we recover the basic reproduction number R₀ by factoring in the history of applied control measures and hence R_0est = R_t/c_est.

The quality of our estimates in Eqs.(10),(11) and (12) depends on the length of the estimation window t_est. Short windows are more sensitive to stochastic fluctuations in incidence, while long windows over-smooth estimates and delay projections [52, 66]. We apply t_est = 5 days, which appears to be a good compromise between the two extremes.