This paper was converted on www.awesomepapers.org from LaTeX by an anonymous user.
Want to know more? Visit the Converter page.

Can the COVID-19 epidemic be controlled
on the basis of daily test reports?

Francesco Casella1 1Francesco Casella is with Dipartimento di Elettronica, Informazione e Bioingegneria, Politecnico di Milano, Piazza Leonardo da Vinci 32, 20133 Milano, Italy; tel. +39 02 2399 3465 francesco.casella@polimi.it
Abstract

This paper studies if and to which extent COVID-19 epidemics can be controlled by authorities taking decisions on public health measures on the basis of daily reports of swab test results, active cases and total cases. A suitably simplified process model is derived to support the controllability analysis, highlighting the presence of very significant time delay; the model is validated with data from several outbreaks. The analysis shows that suppression strategies can be effective if strong enough and enacted early on. It also shows how mitigation strategies can fail because of the combination of delay, unstable dynamics, and uncertainty in the feedback loop; approximate conditions based on the theory of limitation of linear control are given for feedback control to be feasible.

Index Terms:
List of keywords (Control applications, Delay systems, Emerging control applications, Healthcare and medical systems, Modeling)

I Introduction

The first outbreak of the COVID-19 [1] virus epidemic took place in China, starting at the end of 2019, and has since then caused a global pandemic with disruptive effects on public health, social life, and the economy. The pandemic will likely spark a large number of studies to understand its behaviour and to determine effective control strategies.

A wide range of mathematical models have been proposed to describe the dynamic evolution of epidemics, starting from the seminal paper [2], and including a wide range of possibly quite sophisticated models, see e.g. [3] for a comprehensive review. The analysis of these models allows to predict the evolution of the disease over time, its asymptotic behaviour (e.g. endemic disease equilibria vs. eradication), and most importantly how it depends on the model parameters.

Epidemiological models are widely used to design vaccination and treatment strategies based on optimal control, see, e.g. [4] and references therein. They can also be used to design feedback vaccination strategies [5], or even feedback strategies combining different actions such as vaccination, treatment, and culling [6]. Some studies take into the account the feedback effects of behavioural changes in the evolution of an epidemic, see [7] and reference therein.

Most models employed to study control strategies are formulated in terms of ordinary differential equations, e.g. the classical SIR and SEIR models and their variants. In some cases, time delays are also included in the model, to account for the incubation time, see, e.g., [8], [9].

Detailed models of the COVID-19 outbreak have started to appear in the literature. With reference to the outbreak in Italy, [10] proposes an extension of the classic SIR model with eight state variables, while [11] presents a spatially resolved model with nine state variables for each of the 107 provinces of the country. Both models confirmed the appropriateness of the public measures taken by the Italian authorities to contain the virus outbreak. A highly detailed epidemiological model of the UK was used in [12] to predict possible outcomes of the virus outbreak and to suggest the adoption of a suppression policy.

The report [13] attempts to estimate the effects of non-pharmaceutical interventions (NPIs, i.e., public health measures) onto the relative reduction of the reproduction number t\mathcal{R}_{t} of COVID-19, by applying Bayesian methods to data from 11 European countries. Given the estimates of the initial reproduction number R0R_{0}, a reduction by at least 60-70% or more is necessary to suppress exponential growth. The main result of [13] is that lockdown leads to an average reduction of t\mathcal{R}_{t} by 50%, school closure by 20%, other measures around 10%. However, 95% confidence intervals on the reduction factors are huge, e.g. 10% to 80% reduction for lockdown, 0% to 45% reduction for school closure, severely undermining their predictive power. This problem is inherent to the requirement of a large enough data set to be statistically significant, which requires to put countries with very different social habits and very different interpretations of the same measure (e.g. lockdown) in the same data set.

The use of feedback control theory has been advocated early on as a powerful tool to support the management of the COVID-19 outbreak [14]. Unfortunately, most of the existing literature on the control of epidemics involves vaccines or treatments, which are currently not available for COVID-19. Some innovative feedback control strategies have been proposed in preprints at the time of this writing, e.g. [15], which proposes a feedback mitigation strategy based on fast lockdown cycles controlled by a supervisory loop, or [16], advocating a strategy based on massive random testing.

The aim of this paper is to assess the controllability of the COVID-19 outbreak, assuming that the population is sufficiently well mixed and that the decisions of public health measures by the authorities are based on daily reports of positive swab tests, active cases, and total cases. To this aim, a suitably simplified model is presented, which is specifically aimed at capturing the fundamental dynamics of the process which is relevant for feedback control, which turns out to be heavily affected by time delay.

The main result of the analysis is twofold. On one hand, suppression strategies can be effective if enacted early on and with strong enough measures. On the other hand, mitigation strategies turn out to be infeasible if the reproduction number is significantly higher than one, and are in any case limited by the time delays in the feedback loop.

The paper is structured as follows. In Section II, a control-oriented model of the epidemic is introduced and validated against data from the outbreaks in different countries. In Section III, the two above-mentioned strategies are analysed in terms of feedback control. Section IV draws conclusions from the control-theoretical analysis with some recommendations for decision makers and future research.

II Modelling

II-A Derivation of the model

Models of the COVID-19 epidemic such as those mentioned in Sect. I are based on first principles, in the sense that their equations describe the time evolution of different categories of subjects, based on the known mechanisms of infection, recovery, and care of patients. However, their behaviour is ultimately decided by the values of coefficients that must be identified from experimental data, which is preciously scarce in the case of a new disease such as COVID-19. The quality and homogeneity of data used to tune those models are also often questionable: different countries adopt different standards for swab testing, possibly changing them over time; some data get lost because of clerical errors; some countries or regions may report lower numbers than real because of political pressures. Even bona-fide reports may fail to provide reliable data, as revealed by the mismatch between official COVID-19 deaths and additional numbers of deaths on municipal records in previous years. The actual effects of NPIs are still quite uncertain, see [13].

Public policies based on such models cannot thus be applied blindly, but must be adapted and corrected based on the observed outcome. Indicators used by public decision makers include daily reports of a) new positive swab tests, b) current number of infected subjects, and c) total number of reported cases. The crucial question is then: is feedback control feasible at all in such a system?

In order to answer this question, a suitably simplified model of the epidemic is derived here to capture the fundamental dynamics that is relevant for the design of the feedback policy, in particular the dynamic relationship between NPIs and the response of the three above-mentioned indicators. The starting point is the basic SEIR model [3] with the addition of a further compartment LL:

dSdt=βISN\displaystyle\frac{dS}{dt}=-\frac{\beta IS}{N} (1)
dEdt=βISNϵE\displaystyle\frac{dE}{dt}=\frac{\beta IS}{N}-\epsilon E (2)
dIdt=ϵEγI\displaystyle\frac{dI}{dt}=\epsilon E-\gamma{I} (3)
dLdt=γIδL\displaystyle\frac{dL}{dt}=\gamma I-\delta L (4)
dRdt=δL\displaystyle\frac{dR}{dt}=\delta L (5)

were NN is the total population, SS is the number of Susceptible individuals, EE is the number of Exposed individuals, that have caught the infection but are not yet infectious, II is the number of Infectious individuals, LL is the number of subjects which are still iLl, but are no longer infectious due to hospitalization, quarantine, or just because infected subjects are mostly infectious during the first few days after the end of the latency period [17], and RR the number of recovered resistant subjects.

The parameter β\beta accounts for the likelihood of infection per unit time; ϵ\epsilon is the inverse of the average latency time before one becomes infectious, γ\gamma is the inverse of the average time infectious subjects spend by actually infecting other people, and δ\delta is the inverse of the average time subjects remain ill but without infecting others. Given the short time spans involved and the relatively low mortality rate, deaths and births can be neglected, as well as immigration and emigration, that are restricted during the outbreak.

The features of the COVID-19 virus, coupled with the unavailability of effective treatments at the time of this writing, are such that allowing more than a few percent of the population to be infected at any point in time is unacceptable, as doing so would lead to a collapse of the public health system, particularly with reference to the significant fraction of infected subjects needing intensive care to survive the acute respiratory syndrome that the virus can cause. This fact, coupled with the fairly long recovery time (about one month), means that even the worst outbreaks in Western Europe are currently estimated to have infected less than 10% of the population after a few months in the course of the epidemic. This allows to consider S(t)/N1S(t)/N\approx 1 and get rid of Eq. (1). In fact, a significant portion of the population may be not susceptible a priori, e.g. due to genetic reasons. However, absent any concrete evidence of this fact, the precautionary principle suggest to consider the worst case S(t)/N=1S(t)/N=1.

Assuming then a constant value of β\beta, the three eigenvalues of system (2)-(4) are δ-\delta and the two roots pp and rr of

s2+(ϵ+γ)sϵ(βγ)=(sp)(sr).s^{2}+(\epsilon+\gamma)s-\epsilon(\beta-\gamma)=(s-p)(s-r). (6)

If β>γ\beta>\gamma, there is one negative eigenvalue pp and one positive eigenvalue rr. Assuming that the negative exponential mode has already died out, the solution of (2)-(3) is then:

I(t)I(0)ert,\displaystyle I(t)\approx I(0)e^{rt}, (7)
E(t)βI(0)r+ϵert\displaystyle E(t)\approx\frac{\beta I(0)}{r+\epsilon}e^{rt} (8)

The doubling time of I(t)I(t) is Td=log(2)/rT_{d}=log(2)/r. Under the assumption that SNS\approx N, we can approximate the current reproduction number Rt=β/γR_{t}=\beta/\gamma. This formula is not exact during transients, when the number of infectious subjects changes over time, but allows to later obtain some interesting synthetic results, though with some approximation.

In contexts where massive testing was carried out, it was found that about 40% of positive tested subjects are entirely asymptomatic, despite being infectious [18], [19], which makes COVID-19 particularly insidious. This suggests that a similar fraction of the infectious subjects goes unnoticed and escapes the testing process in the general population. Furthermore, swab tests are affected by more than 20% false negatives [20], and some mildly symptomatic subjects may also end up not being tested. The ratio α\alpha between the number the infectious subjects I(t)I(t) at a certain time tt and the number of infectious subjects It(t)I_{t}(t) at the same time tt that will eventually get tested positive is then likely between two and three. On the other hand, α\alpha is only relevant to determine when the ratio S(t)/NS(t)/N starts decreasing significantly below one, providing some degree of herd immunity, as all other important indicators, namely the mortality ratio, the ratio of hospitalized subjects and the ratio of subjects requiring intensive care, are all referred to It(t)I_{t}(t).

Assuming that α\alpha is constant, one can use Eqs. (2)-(5) to also describe the dynamics of the fraction of subjects that are eventually tested positive through the various stages of the disease, Et(t)E_{t}(t), It(t)I_{t}(t), and Lt(t)L_{t}(t).

In most cases, subjects are only tested after they show serious symptoms, which happens on average τt\tau_{t} days after they have become infectious. The lab processing also introduces a delay τr\tau_{r} before reports are available. Although in principle it is possible to provide the results of the test in a few hours, the average reporting time is usually much longer because of limited equipment availability, up to one week or more.

The NPIs mentioned earlier (lockdown, school closures, etc.) reduce the rate of infection β\beta, hence the current reproduction number tβ/γ\mathcal{R}_{t}\approx\beta/\gamma. These measures are varied and can be applied progressively. We can then assume that the time-varying parameter β\beta is in fact a function of a representative manipulated variable u(t)u(t), with uu indicating the intensity of adopted public health measures on a scale from 0 (no intervention) to 1 (full lockdown and isolation of all individuals). The β()\beta(\cdot) function is thus monotonously decreasing from the value β0\beta_{0}, when no social restrictions are enforced, to zero, corresponding to the total isolation of each person in the contry. Of course β\beta is also a function of other unknown factors d(t)d(t) that act as disturbances on the system, e.g. mutations of the virus or changes in social behaviour which are not directly mandated by the authorities. Considerable uncertainty is involved in the estimation of the effects of different interventions in terms of reduction of β\beta or, equivalently, of t\mathcal{R}_{t}, see [13], hence β()\beta(\cdot) is also uncertain.

The control-oriented model can thus be formulated as a state-space system with output delays:

dEt(t)dt\displaystyle\frac{dE_{t}(t)}{dt} =β(u(t),d(t))It(t)ϵEt(t)\displaystyle=\beta(u(t),d(t))I_{t}(t)-\epsilon E_{t}(t) (9)
dIt(t)dt\displaystyle\frac{dI_{t}(t)}{dt} =ϵEt(t)γIt(t)\displaystyle=\epsilon E_{t}(t)-\gamma I_{t}(t) (10)
dLt(t)dt\displaystyle\frac{dL_{t}(t)}{dt} =γIt(t)δLt(t)\displaystyle=\gamma I_{t}(t)-\delta L_{t}(t) (11)
dTt(t)dt\displaystyle\frac{dT_{t}(t)}{dt} =ϵEt(t)\displaystyle=\epsilon E_{t}(t) (12)
Nr(t)\displaystyle N_{r}(t) =ϵEt(tτm)\displaystyle=\epsilon E_{t}(t-\tau_{m}) (13)
Ar(t)\displaystyle A_{r}(t) =It(tτm)+Lt(tτm)\displaystyle=I_{t}(t-\tau_{m})+L_{t}(t-\tau_{m}) (14)
Tr(t)\displaystyle T_{r}(t) =Tt(tτm),\displaystyle=T_{t}(t-\tau_{m}), (15)

where β(u,d)\beta(u,d) is an uncertain function, ϵ\epsilon, γ\gamma, δ\delta are uncertain constant parameters, τt\tau_{t}, τr\tau_{r} are uncertain parameters, τm=τt+τr\tau_{m}=\tau_{t}+\tau_{r} is the overall measurement delay, and Tt(t)T_{t}(t) computes the cumulative number of positive tested subjects.

The first measured variable of the process is the number of new daily reported cases Nr(t)N_{r}(t), which is affected by the overall delay τm\tau_{m}. The second measured variable is the number of reported active cases Ar(t)A_{r}(t), i.e. the number of subjects for which a positive test report has been received and a double negative test has not yet been issued to certify their recovery. As τt\tau_{t} and τr\tau_{r} are similar (several days), 2τrτt+τr=τm2\tau_{r}\approx\tau_{t}+\tau_{r}=\tau_{m}, leading to Ar(t)=It(tτm)+Lt(tτm)A_{r}(t)=I_{t}(t-\tau_{m})+L_{t}(t-\tau_{m}). The third measured variable is the total cumulated number of reported positive swab test reports Tr(t)=Tt(tτm)T_{r}(t)=T_{t}(t-\tau_{m}).

Note that the model (10)-(15) has a time delay on the output equations. Since input/output dynamics only will be considered in the next Section, an equivalent representation could be used where the delay is applied to the input instead.

II-B Validation and Tuning

The goal of the model is to describe the dynamic response of the Nr(t)N_{r}(t), Ar(t)A_{r}(t), and Tr(t)T_{r}(t) indicators to the application of NPIs by central authorities, described by changes in u(t)u(t). Four outbreaks cases were selected, all characterized by step changes of u(t)u(t) at the central government level, in order to make the validation easier: China, with data taken from [1], Italy, France, and the UK, with data taken from [21], which reports data from national authorities.

In the case of Italy and UK, some partial restrictions were introduced first, causing a noticeable delayed reduction of the exponential increase rate of new cases, then a full lockdown was prescribed, whose effect was to change the positive exponential growth of new cases into a negative exponential decay after some delay. We then assume that β=β0\beta=\beta_{0} at t=0t=0, then β\beta undergoes a step reduction to β=ρ1β0\beta=\rho_{1}\beta_{0} at t=t1t=t_{1} and a further step reduction to β=ρ2β0\beta=\rho_{2}\beta_{0} at time t=t2t=t_{2}. In the case of China and France no NPIs were enforced before full lockdown, yet a reduction of the exponential growth rate well ahead of the effects of lockdown is clearly visible, possibly due to social feedback effects or other disturbances. Two step reduction were thus applied also in those cases.

The parameters of the model were tuned manually to obtain a good fit with the available data. In particular, τm\tau_{m} is easily tuned to match the delay between the lockdown and the peak of Nr(t)N_{r}(t), t=β/γ\mathcal{R}_{t}=\beta/\gamma and ρ1\rho_{1} determine the growth rates in the pre-lockdown behaviour of Nr(t)N_{r}(t), ρ2\rho_{2}, γ\gamma and ϵ\epsilon determine the shape and decay rate of the post-lockdown behaviour of Nr(t)N_{r}(t). Tr(t)T_{r}(t) is the integral of Nr(t)N_{r}(t), so fitting it helps refining the parameter tuning considering the noisy nature of Nr(t)N_{r}(t), which is also affected by a weekly oscillation due to repeating lab schedules. Finally, δ\delta is tuned to match the peak of active cases Ar(t)A_{r}(t), which is much wider and delayed than the peak of Ir(t)I_{r}(t). A detailed error analysis was not performed and could be the subject of future work considering a more extensive data set; in any case, high parameter accuracy is not required to support the forthcoming analysis.

The values obtained for the four outbreaks are reported in Table I. They are fairly consistent with each other and are compatible with the ones reported in [11] and [21]. The initial reproduction number 0\mathcal{R}_{0} and the current reproduction numbers 1\mathcal{R}_{1} and 2\mathcal{R}_{2} computed after each change of β\beta are reported, as well as the doubling times Td0T_{d0} and Td1T_{d1} of the unstable mode before and after the first change of β\beta, and the maximum reproduction number l\mathcal{R}_{l} corresponding to the controllability limit derived in Section III-B4.

TABLE I: Model parameters (time constants in days)
Period β0\beta_{0} 1/γ1/\gamma 1/ϵ1/\epsilon 1/δ1/\delta t1,ρ1t_{1},\rho_{1} t2,ρ2t_{2},\rho_{2} 0\mathcal{R}_{0} 1\mathcal{R}_{1} 2\mathcal{R}_{2} Td0T_{d0} Td1T_{d1} τm\tau_{m}
China 18/01/2020 – 11/02/2020 1.6 2.5 5.0 N/A 0, 0.63 5, 0.160 4.0 2.5 0.64 2.5 4.3 12
Italy 22/02/2020 – 01/05/2020 1.3 3.1 4.3 33 2, 0.56 19, 0.205 4.0 2.5 0.82 2.6 5.2 9
France 28/02/2020 – 03/05/2020 1.3 2.9 5.0 29 2, 0.60 17, 0.195 3.8 2.28 0.74 2.8 5.3 12
UK 01/03/2020 – 10/05/2020 1.28 2.8 6.2 N/A 15, 0.65 25, 0.270 3.6 2.34 0.97 3.4 5.8 10
Refer to caption
Figure 1: Italian outbreak validation: new daily reported cases
Refer to caption
Figure 2: Italian outbreak validation: active and total reported cases

The detailed results for the case of Italy, computed with the code available online in [22], are reported in Figs. 1-2. The interplay between EtE_{t} and RtR_{t} accounts for the growth or decay of cases, depending on t\mathcal{R}_{t}. The EE compartment is essential to explain why Nr(t)N_{r}(t) does not drop sharply, when β\beta is sharply reduced, once the τm\tau_{m} delay has elapsed, while the LL compartment is necessary to explain the much slower dynamic response of active reported cases Ar(t)A_{r}(t). The validation results of the other cases are reported in [23].

III Control

The effects of the application of the two control policies outlined in the Introduction will now be analysed. The title of this section may well be ”Respect the Unstable” [24]: feedback control strategies should not be applied light-heartedly to safety-critical unstable systems.

III-A Suppression

The suppression strategy can be brutally summarized in the following terms: as soon as Ar(t)A_{r}(t) reaches a value AsA_{s} which is scary enough to decision makers to overcome their reluctance to disrupt the social and economic life of their country, drastic containment measures are taken:

u(t)={0,Ar(t)<Asu¯,Ar(t)As.u(t)=\left\{\quad\begin{matrix}0,&A_{r}(t)<A_{s}\\ \bar{u},&A_{r}(t)\geq A_{s}\end{matrix}\right.. (16)

If the threshold AsA_{s} is crossed at time tst_{s} and u¯\bar{u} is large enough, then β(u¯)/γ<1\beta(\bar{u})/\gamma<1 and thus r<0r<0. Assuming also d(t)d(t) remains constant for ttst\geq t_{s}, Eqs. (9)–(11) form a homogeneous LTI system with three negative eigenvalues rr, pp, and δ-\delta. The actual number of eventually tested positive exposed subject Et(t)E_{t}(t) will start decaying immediately; however, the number of new daily reported infectious cases Nr(t)N_{r}(t) will continue its exponential growth for τm\tau_{m} days, before starting to decay as well. The number of reported active cases Ar(t)A_{r}(t), hence the required number of beds in hospitals and intensive care units, will also stop increasing exponentially after τm\tau_{m} days, but will continue growing and peak much later, due to the much slower time constant 1/δ1/\delta, see, e.g., Fig. 2. Then, states and outputs asymptotically approach the equilibrium in the origin, that corresponds to the eradication of the virus.

The peak value of active cases Ap=maxAr(t)A_{p}=\max A_{r}(t) can be computed by numerical integration of Eqs. (9)–(11) and (14). The ratio M=Ap/Ar(ts)M=A_{p}/A_{r}(t_{s}) can be quite large, e.g. M=8M=8 for the Italian outbreak and M=10M=10 for the French outbreak.

Assuming that a fraction σ\sigma (about 4%4\% in Italy) of active cases requires intensive care, and that NicN_{ic} intensive care beds are available, a wise choice of AsA_{s} requires σAr(t)<Nict\sigma A_{r}(t)<N_{ic}\enspace\forall t; hence, As<Nic/σMA_{s}<N_{ic}/\sigma M. Political decision-makers without a training in mathematical modelling may have difficulties in understanding the role and magnitude of factor MM and may be caught by surprise once it is too late to act.

III-B Mitigation

III-B1 Policy statement

The basic idea of mitigation policies is to manage the outbreak, in particular the trajectory of Ar(t)A_{r}(t), in order to avoid overloading the public health system, without trying to suppress it. This strategy was followed until 16 March 2020 by the UK government, which aimed at achieving herd immunity [25], and until at least 10 June 2020 by the Swedish government [26].

III-B2 Mathematical formalization

The first step to enact this strategy is to compute a reference control policy u0(t)u^{0}(t), obtained by the application of suitable NPIs over time, whose effects on β\beta is accurately calibrated, leading to reference trajectories Nr0(t)N^{0}_{r}(t), Ar0(t)A^{0}_{r}(t) and Tr0(t)T^{0}_{r}(t) for the corresponding indicators, respecting the constraint σAr0(t)<Nic\sigma A_{r}^{0}(t)<N_{ic}. These trajectories can be obtained by means of constrained optimal control, using sophisticated models of the epidemic as, e.g., the ones reported in [12], [11], and [10].

The unstable nature of the state trajectories while r>0r>0 makes an open-loop implementation of this policy infeasible, unless one wants to risk runaways that can cause the collapse of the public health system. The reference trajectory should rather be followed by adapting the public policy measures u(t)u(t) in real time, based on the values of Nr(t)N^{r}(t), Ar(t)A^{r}(t), and Tr(t)T^{r}(t), which are constantly monitored by the authorities. This corresponds in principle to closing a feedback loop to stabilize the unstable reference trajectory, see Fig. 3.

Refer to caption
Figure 3: Mitigation control architecture

III-B3 Ideal Feedback controller design

The process model (9)-(14) can be linearized around the reference trajectory at time tat_{a}, obtaining a linear model with constant coefficients except for the terms β(u0(ta),d(ta))\beta(u^{0}(t_{a}),d(t_{a})) and β(u0(ta),d(ta))u\frac{\partial\beta(u^{0}(t_{a}),d(t_{a}))}{\partial u}, which depend on tat_{a} for non-trivial reference control trajectories u0(t)u^{0}(t). For the sake of the subsequent analysis, we assume that these parameters change over a time scale which is much longer than the time scale of the closed-loop system feedback response, a common assumption when dealing with gain-scheduling control, and thus consider them as constants, with the value they have at time tat_{a} around which the feedback stability analysis is performed. The transfer functions of the linearized process then reads:

ΔNr(s)Δu(t)\displaystyle\frac{\Delta N_{r}(s)}{\Delta u(t)} =μ(ta)(s+γ)(sp(ta))(sr(ta))eτms\displaystyle=\frac{\mu(t_{a})(s+\gamma)}{(s-p(t_{a}))(s-r(t_{a}))}e^{-\tau_{m}s} (17)
ΔAr(s)Δu(t)\displaystyle\frac{\Delta A_{r}(s)}{\Delta u(t)} =μ(ta)(s+(γ+δ))(sp(ta))(sr(ta))(s+δ)eτms\displaystyle=\frac{\mu(t_{a})(s+(\gamma+\delta))}{(s-p(t_{a}))(s-r(t_{a}))(s+\delta)}e^{-\tau_{m}s} (18)
ΔTr(s)Δu(t)\displaystyle\frac{\Delta T_{r}(s)}{\Delta u(t)} =μ(ta)(s+γ)(sp(ta))(sr(ta))seτms\displaystyle=\frac{\mu(t_{a})(s+\gamma)}{(s-p(t_{a}))(s-r(t_{a}))s}e^{-\tau_{m}s} (19)

where μ(ta)=ϵβ(u0(ta),d(ta))uI0(ta)\mu(t_{a})=\epsilon\frac{\partial\beta(u^{0}(t_{a}),d(t_{a}))}{\partial u}I^{0}(t_{a}), and p(ta)p(t_{a}) and r(ta)r(t_{a}) are the eigenvalues of system (9)-(10) linearized at t=tat=t_{a} around the reference trajectory.

By making the very optimistic assumptions that the parameters ϵ\epsilon, γ\gamma, δ\delta, and τm\tau_{m} are constant and perfectly known, that the function β(u,d)\beta(u,d) that expresses the effects of public policy decisions is perfectly known, monotonously decreasing and smooth with respect to uu, and that d(t)0d(t)\equiv 0, one can design the feedback controller CFB\mathrm{C_{FB}} as a linear controller C(s)C(s) with gain scheduling, that compensates for the nonlinearity of the process gain, resulting in a linear and (approximately) time-invariant loop dynamics. While doing so, one should also account for an additional delay τc\tau_{c} of 2÷42\div 4 days within the controller, corresponding to the decision making and implementation delay. Assuming one wants to use Nr(t)N_{r}(t) for feedback control, the overall control law is thus:

u(t)=u0(t)+1μ(tτc)uf(tτc)\displaystyle u(t)=u^{0}(t)+\frac{1}{\mu(t-\tau_{c})}u_{f}(t-\tau_{c}) (20)
uf(s)=C(s)[Nr0(s)Nr(s)]\displaystyle u_{f}(s)=C(s)\left[N_{r}^{0}(s)-N_{r}(s)\right] (21)

and the loop transfer function of the controlled system is

L(s)=C(s)μr(s+γ)(sp(tτc))(sr(tτc))es(τm+τc),L(s)=C(s)\frac{\mu_{r}(s+\gamma)}{(s-p(t-\tau_{c}))(s-r(t-\tau_{c}))}e^{-s(\tau_{m}+\tau_{c})}, (22)

where μr\mu_{r} is the ratio between the actual value of the gain μ\mu of transfer function (17) and its reference value used for gain scheduling. In ideal conditions, μr=1\mu_{r}=1, though results from [13] imply this gain is subject to significant uncertainty. In case one wants to control Ar(t)A_{r}(t) or Tr(t)T_{r}(t), similar consideration apply, using (18) or (19) instead of (17).

In all the three cases, the loop transfer function reveals the very dangerous nature of this process, which features a time delay τ\tau, an uncertain gain μr\mu_{r}, and an unstable pole with time constant TT if β(ta)/γ>1\beta(t_{a})/\gamma>1, where

T=1r=Tdlog(2)\displaystyle T=\frac{1}{r}=\frac{T_{d}}{log(2)} (23)
τ=τt+τr+τc.\displaystyle\tau=\tau_{t}+\tau_{r}+\tau_{c}. (24)

Considering the values reported in Table I, Td=4.3÷5.8T_{d}=4.3\div 5.8 before lockdown, while τ=9÷12\tau=9\div 12 days.

III-B4 Control feasibility

In order to guarantee some robustness of the system performance against the large gain uncertainty of the process, the Bode plot of |L(jω)|\mathinner{\!\left\lvert L(j\omega)\right\rvert} should maintain a roughly constant slope over a sufficiently wide interval around the crossover frequency ωc\omega_{c}, thus approximating Bode’s ideal loop transfer function.

When t(t)>1\mathcal{R}_{t}(t)>1, hence r(t)>0r(t)>0, the analysis reported in [27], Sect. 4.6, leads to conclude that if one wants to limit the maximum norm of the sensitivity function Ms<2M_{s}<2 to achieve some robustness, the product of the unstable pole rr and of the time delay τ\tau should be rτ<0.156r\tau<0.156. Introducing the doubling time TdT_{d}, this condition becomes:

τTd<0.225,\frac{\tau}{T_{d}}<0.225, (25)

i.e., under very optimistic assumptions on the knowledge of the process parameters, feedback control is feasible only if the overall loop delay is less than a quarter of the doubling time of the outbreak. Considering the values reported in Table I, this limitation translates into t(t)<1.1\mathcal{R}_{t}(t)<1.1.

When t(t)<1\mathcal{R}_{t}(t)<1, hence r(t)<0r(t)<0, [27] concludes that the maximum crossover frequency is 1.57/τ1.57/\tau; hence, the time constant of the response of the feedback controller to unexpected disturbances cannot be less than T=τ/1.57T=\tau/1.57.

Unfortunately there is no theorem that can be directly invoked to prove that any feedback control policy would not suffer from the same limitations of a carefully scheduled linear controller. However, the results of this analysis provides an insightful benchmark, pointing out the crucial role of measurement and decision delays, which should explicitly be taken into account for feedback control design, and minimized as much as possible. It also suggest that robust feedback control may not be feasible around trajectories where t\mathcal{R}_{t} is significantly above one. For example, this may explain the runaway scenarios happening with non-negligible likelihood in [15], when taking into account the statistical distribution of the uncertain process parameters.

IV Conclusions

Governments all the world over are faced with very challenging life-or-death decisions regarding the management of the COVID-19 epidemic, involving the balance between public health and economic issues. In order to take such decisions, they rely on expert advice based on the results of epidemiological mathematical models and on daily case reports, based on swab test results.

This paper puts the problem in a control systems perspective as a feedback control problem, using a simple model to capture the control-relevant dynamic response of those reports to the application of NPIs. The model was tuned and validated with data from four different outbreaks.

These are the main results of the analysis:

  • The suppression strategy is effective if NPIs are strong enough to obtain Rt<1R_{t}<1, but it requires to understand the role of the multiplicative factor MM to correctly decide when it is the right time to enforce them.

  • Mitigation strategies are limited by the combination of delay, uncertainty, and unstable dynamics. Designing robust stabilizing controllers around trajectories with t>1.1\mathcal{R}_{t}>1.1 is likely to be difficult or impossible. Reducing the overall delay by 50% would bring the limit to t>1.2\mathcal{R}_{t}>1.2. This information is particularly relevant for the management of the reopening phase after lockdown.

  • Measurement and decision delays play a crucial role in determining the feedback control performance and stability; hence, they should be explicitly taken into account in the design of any feedback controller, and minimized as much as possible, e.g., by promoting fast testing policies and technologies.

  • The analysis and design of NPIs can benefit from control theory tools, possibly suggesting viable solutions or pointing out shortcomings of proposed strategies, that are not obvious to epidemiologists and physicians.

At the time of this writing, the emerging consensus seems to be that the safest policy to address exponentially growing COVID-19 outbreaks is to apply aggressive enough suppression policies; the results reported in this paper can further motivate why this is actually the case.

These results could also be useful to devise effective and safe strategies to cope with the reopening phase that countries face after successfully suppressing the first outbreak, in particular if the number of new daily cases becomes too large to allow for testing, tracing and tracking of individual cases.

References

Supplementary material

Detailed results of the model validation in the cases of China, France, and UK are discussed in this appendix, which does not fit the strict 6-page limit of IEEE Control Systems Letters publications, and therefore is only made available in the arXiv version of the paper.

IV-A China

Fig. 4 reports the validation of the Nr(t)N_{r}(t) trajectory for the Chinese outbreak, based on data from [1].

The lockdown became effective on day 5 for the city of Wuhan, which had the majority of cases at the time, and on day 6 for all the Hubei region. Since it is not possible from available data to figure out what are the separate effects of two such very close events, the second one was lumped together with the first one, which produced the largest number of reported cases.

The data at the beginning of the time series until day 5 doesn’t really fit well an exponential curve, but that is probably due to issues with the collection of swab test results. However, the effect of the time delay τm\tau_{m} is very clearly visible, as the number of new reported daily cases only stopped growing exponentially after 12 days. The behaviour of Nr(t)N_{r}(t) after the peak allowed to tune the γ\gamma and ϵ\epsilon parameters, and to estimate a reproduction number t=0.64\mathcal{R}_{t}=0.64 immediately after the effect of lockdown was felt.

Refer to caption
Figure 4: Chinese outbreak validation: new daily reported cases

IV-B France

Fig. 5 reports the validation of the Nr(t)N_{r}(t) trajectory for the French outbreak, based on data from [21].

As already mentioned in Sect. II-B, the doubling time increased significantly after day 14, corresponding to a change of β\beta happened around day 2, given the delay τ=12\tau=12. The cause of this change is unknown, it may actually have been caused by some feedback effects on social behaviour caused by the growing fear of contagion sparked by the nearby Italian case. The trend of Nr(t)N_{r}(t) after the delayed effect of the lockdown on day 17 is very noisy but clearly decreasing on average. This is confermed by looking at the very good match in the validation of the Tr(t)T_{r}(t) trajectory reported in Fig. 6: as Tr(t)T_{r}(t) is the integral of Nr(t)N_{r}(t) (see (15)), high-frequency noise is filtered out.

As to the trajectory of active cases Ar(t)A_{r}(t) shown in Fig. 6, the matching of the simulation results with the data is very good for the first 25 days, then some mechanism seems to have somewhat dampened the experimental curve more efficiently than the model in this paper. In any case, the model overestimates the peak value by around 15%, which is not that bad, considering how simple the model is.

Refer to caption
Figure 5: French outbreak validation: new daily reported cases
Refer to caption
Figure 6: French outbreak validation: active and total reported cases

IV-C UK

Finally, we take into consideration the case of UK, again using data taken from [21]. The virus was basically allowed to run unchecked until March 16, 2020 (day 15), when the government issued some recommendations to avoid going to pubs or theaters, without however closing them or issuing any legally binding rules or restriction. This seems to have slowed down the exponential growth after a 10-day delay, increasing the doubling time from 3.4 to 5.8 days.

Refer to caption
Figure 7: UK outbreak validation: new daily reported cases

Full lockdown was announced on Mar 23 and went into effect on Mar 26 (day 25). Also in this case, the effect on Nr(t)N_{r}(t) was felt after a significant delay, in this case of about 10 days, see Fig. 7. Once that time interval was elapsed, the number of new daily positive swab tests plateaued, instead of starting to decay, settling down on a trajectory with estimated t=0.97\mathcal{R}_{t}=0.97, very close to the stability limit.

Also in this case, the experimental data about Nr(t)N_{r}(t) are very noisy; however, an excellent match with data is obtained with the cumulated Tr(t)T_{r}(t) trajectory, which is less sensitve to nois because of the additional integral effect, see Fig. 8. Unfortunately, data for Ar(t)A_{r}(t) were not available from [21], so it was not possible to validate the active cases output in this case.

Refer to caption
Figure 8: UK outbreak validation: active and total reported cases