key: cord-0529553-hi4ifty0
authors: Burke, Kevin; Barmish, B. Ross
title: A Data-Driven Control-Theoretic Paradigm for Pandemic Mitigation with Application to Covid-19
date: 2020-08-14
journal: nan
DOI: nan
sha: b1b5053ad5a6cd2d1a8b6e6a426290d8accb31bf
doc_id: 529553
cord_uid: hi4ifty0

In this paper, we introduce a new control-theoretic paradigm for mitigating the spread of a virus. To this end, our discrete-time controller, aims to reduce the number of new daily deaths, and consequently, the cumulative number of deaths. In contrast to much of the existing literature, we do not rely on a potentially complex virus transmission model whose equations must be customized to the"particulars"of the pandemic at hand. For new viruses such as Covid-19, the epidemiology driving the modelling process may not be well known and model estimation with limited data may be unreliable. With this motivation in mind, the new paradigm described here is data-driven and, to a large extent, we avoid modelling difficulties by concentrating on just two key quantities which are common to pandemics: the doubling time, denoted by $d(k)$ and the peak day denoted by $theta(k)$. Our numerical studies to date suggest that our appealingly simple model can provide a reasonable fit to real data. Given that time is of the essence during the ongoing global health crisis, the intent of this paper is to introduce this new paradigm to control practitioners and describe a number of new research directions suggested by our current results.

The main objective in this paper is to introduce a new control-theoretic paradigm for mitigation of the spread of pandemic disease. In the sequel, the discrete-time control variable u(k) corresponds to the "degree of mitigation" which could be clinical or nonclinical in nature. The motivation for the paradigm to follow is simple to explain: When a new pandemic emerges on the world scene, standard virus transmission models are typically adopted in an effort to predict and control the uncertain future. However, with limited epidemiological information available for a novel virus, existing models may be inappropriate; i.e., such a virus may have many features in its transmission and mitigationresponse dynamics which are not well understood and not captured using virus models from the past. Accordingly, with this motivation in mind, data-driven modelling is the focal point in this paper.

In the existing literature, the majority of papers on virus transmission are based on the so-called "SIR" compartmental model [1] which, at a given point in time, places individuals into one of three disjoint classes: the susceptible S class, the infectious I class, and the recovered R class. The wide variety of extensions to this model include considerations such as death due to the disease, the latency time for infected individuals to become infectious, severity of illness, quarantining individuals, spatial and demographic effects, and social network structure; e.g., see the detailed accounts in [3, 2, 4] . As stated in [5] , there are "a thousand and one epidemic models," and, therefore, customizing the model to the specifics of the situation at hand is non-trivial. This "customization" issue is epitomized by the variety of models even within the Covid-19 context, for example, in [6, 7, 8, 9] , we see a range of 7-16 possible states emphasizing different epidemiological considerations. In the early days of a pandemic, with limited data, not only is it difficult to build the structural model, but also to estimate its parameters. A particularly thorny modelling issue is that many infected individuals go undetected [10] , in large part due to asymptomatic or mild cases. This runs counter to many SIR modelling efforts which assume that the detected cases are equivalent to the true cases when, in fact, the former will underestimate the latter.

Our approach in this paper, largely intended to counter concerns along the lines above, is motivated by data reliability issues which arise in the literature. To this end, our starting point is to work with data in the form of death numbers, which are typically higher fidelity as these are much less likely to be undetected or under-reported; e.g., see [11] . Although at first sight, the SIR-related models of [12] and [13] also appear only to use death data, in fact, their models borrow infection and death-rate parameters from other papers in the literature. In another recent paper [14] , this issue is partially remedied by using both death and case data with some additional parameters to account for the difference between detected and true cases. In practice, the high variability in the detection process seems consistent with the wide uncertainty bands produced by [13] for the proportion of infected individuals in the population.

In summary, while research efforts aimed at understanding virus transmission mechanisms are valuable and important, the challenges in modelling will not be lost on the reader. Therefore, as mentioned above, we focus entirely on the death process data and propose a parsimonious three-parameter model. We expect that this approach will appeal to the control community wherein epidemiology is not a core competency. That being said, although we avoid transmission dynamics, we make use of two quantities which are common to epidemics over limited time frames: the so-called "doubling time" and "peak day" for deaths; over longer horizons, we note that periodic peaking may occur [15] . Thus, our approach can be viewed as lying between the epidemiologically-based SIR models and empirically-based phenomenological models in [11] and [16] .

In addition to the considerations above, in this paper, we also include the dynamics of a controller reflecting the effect of mitigation measures aimed at reducing new deaths to zero. To the best of our knowledge, the earliest work bringing control-theoretic methodologies to epidemiological modelling dates back to the 1970s; see [17, 18, 19] . Building on this, in [20] the theoretical properties of various mitigation techniques are considered, while others focussed more specifically on vaccination [21] , the combination of vaccination and treatment [22] , non-clinical interventions such as social distancing, quarantining, and education [23] , system delay [24, 25] and spatiotemporal effects [26] ; more recently, the Covid-19 outbreak has been emphasized in [27, 28, 29] . Interestingly, all of these control strategies assume underlying SIR-type models which, as previously discussed, present various challenges in practice. In contrast to these approaches, a key feature of our datadriven approach is the error dynamics involving the comparison of model-predicted deaths to actual deaths. This enables both the model and the mitigation controller to be adapted over time.

The remainder of the paper is organised as follows: In Section 2 we describe the preliminaries involving the use of data and error dynamics at a high level. Subsequently, in Sections 3 and 4, the details of our new model and its qualitative properties are provided for the case of constant mitigation. In Section 5, we describe the process of model parameter estimation and illustrate its use via numerical examples involving Covid-19 data from Brazil and Mexico, two of the countries among those with the highest of death rates as of mid-2020. Finally, in Section 6, conclusions and directions for future research are described with emphasis being on issues of a control-theoretic nature.

In the sequel, we take k to be the index indicating the day number and u(k) to be the corresponding level of mitigation provided by the controller. In this first paper aimed at introducing our new paradigm, we do not consider the detailed mechanics of mitigation. Suffice it to say, large values of u(k) might represent stronger mitigation measures such as government mandates on social distancing and the use of masks and therapeutics whereas smaller values might correspond to "relaxation" of the rules.

Starting at k = 0, at a general level, one begins with an equation for new deaths

where the cumulative total deaths are naturally constrained to be D(k) = D(k−1)+N(k).

We propose an attractive form for f (D, N, u) in Section 3, based on the doubling time and peak day quantities from epidemiology.

Our control-theoretic paradigm begins with the acquisition of daily death data which is collected to obtain estimates of the parameters describing the function f above. Once the model is fixed, consistent with the tenets of receding horizon control, one can make predictions to determine if the control sequence u(k) is mitigating in the sense that lim k→∞ N(k) = 0. In practice, the speed of convergence is also a concern and strictly reaching the zero limit may not be required if the number of deaths is an acceptably small fraction of the population size. As a pandemic unfolds over time, the model estimation should be updated periodically as new data is obtained. That is, letting N a (k) denote the "actual" number of new deaths on day k, as depicted in Figure 1 , we use the error e(k) .

= N a (k)−N(k) to drive the update of the model, the predicted deaths, and the associated adaptation of the controller u(k).

There are many ways one could proceed when using the error; e.g., one can down-weight the distant past or smooth the noise using cumulative errors. However, such considerations are beyond the scope of this article. 

Our analysis begins with the so-called doubling time, a widely reported metric used in practice to characterise the number of days d taken for the new deaths to double. * If, for * Note our use of the doubling time for new deaths whereas some authors and data providers use this terminology for cumulative total deaths. example, the doubling time is constant, the number of new deaths on day k + 1 is given by N(k + 1) = 2 1/d N(k) so that N(k + d) = 2N(k) as expected. In practice, d = d(k) is time-varying for reasons such as immunity building up in the population, changes to mitigation strategies, and medical developments. Therefore, with initial values d(0) = d 0 and N(0) = N 0 , viewed as "baseline" quantities at the point from which we model, we obtain new deaths as

In the equation above d(k) can be quite general in its functional form. We now specialise along the lines described in Section 1. That is, we structure d(k) so that it is consistent with many standard epidemic models which exhibit "peaking" behaviour. Thus, with d 0 > 0, the count on new deaths N(k) climbs until some peak day

which may be time-varying due to changes in the level of mitigation characterized by the controller u(k). After the peak day, d(k) becomes negative so that we see declining N(k); in this post-peak regime, |d(k)| may be viewed as the halving time.

We capture the peaking behaviour by specifying the doubling model as

which increases as k approaches θ(k) > 0 from the left. If, for some k * , we have k < k * < θ(k) and k > k * > θ(k), then a single peak occurs at k = k * = θ(k * ). More generally however, θ(k) may be referred to as the "anticipated" peak day since it can move further or closer in time due to the variations in u(k) which could even yield multiple peaks.

An important starting point in this framework is the case when the level of mitigation is being held fixed over some extended time period. To this end, we consider the controller u(k) = u 0 for all k. In practice, this would correspond to way the model is typically used; i.e, the mitigation level is held fixed between model updates and associated assessments of the efficacy of control measures in place. Subsequently, if the updated model parameters predict a worsening prognosis, decision makers might increase the level of mitigation u(k). For the "unit step" input control above, the corresponding peak day is denoted

and we obtain the equation for new deaths

which is readily solved for k ≥ 0. Indeed, expressing N(k) as a product followed by summation of exponents yields

The appealingly simple solution for N(k) above lends itself to various insights about the evolution and qualitative behavior of deaths. First, it is evident by inspection that N(k) is bell-shaped due to the negative k 2 exponent, and this is clear from Figure 2 . It is also straightforward to study the dependence of N(k) on θ 0 and d 0 by first noting that

For the most important case where d 0 and θ 0 are positive, we see that N(k) increases in θ 0 , decreases in d 0 for k < k N 0 , and increases in d 0 for k > k N 0 .

It is apparent that k N 0 represents the number of days until the pandemic reverts back to the early stage baseline level; i.e., N(k N 0 ) = N 0 . However, in contrast to the early stage, for k > k N 0 , new deaths N(k) are now on a downward trajectory which in some sense can be viewed as signalling the end point of the pandemic. Finally, it is also of interest to characterize the number of deaths on the peak day, N(θ 0 ) = 2 θ 0 +1 2d 0 N 0 which, clearly, increases in θ 0 and decreases in d 0 .

The peak N(θ 0 ) above is important in that it correlates strongly with "anticipated pressure" on hospitals; e.g., it is a predictor of stress on resources such as intensive care units. In particular, given the concern that N(θ 0 ) exceeds some critical level N max , one can easily use the equations above to study safety margins such as

which indicates how far from N max , as a percentage, peak deaths will be. When S < 0, since this corresponds to N(θ 0 ) > N max , one might consider increasing the level of mitigation by updating u(k) to reduce a potential crisis. In Figure 3 , we display a realm of possible outcomes for this safety margin where, for example, with an initial doubling time of 10 days, an anticipated peak of up to about 14 days, the (d 0 , θ 0 ) point lies in the green-colored safety zone. 

The total number of deaths, given by

is not summable as a closed form solution. However, many properties of D(k), inherited from N(k), are nevertheless clear: First, since each N(i) term is a point on a bell-shaped curve, D(k) must be a nondecreasing sigmoidal function. Moreover, since N(i) is increasing in θ 0 for all i, D(k) must also be increasing in θ 0 . This makes the importance of reducing θ 0 by mitigation quite clear. That is, a reduced θ 0 leads to a smaller peak which happens earlier, and, in turn, this lowers the total number of deaths. Furthermore, it is also apparent that D(k) decreases with respect to d 0 for k ≤ k N 0 . For k > k N 0 , the N(i) terms are increasing in d 0 which makes the global behaviour of D(k) with respect to d 0 non-trivial. However, recalling that k N 0 can be viewed as signalling the end of the pandemic, for reasonable d 0 , the terms entering D(k) beyond N(k N 0 ) will be quite small; i.e., the portion of D(k) which increases with respect to d 0 will be small enough such that D(k) decreases with d 0 . That being said, for very large d 0 , it is easy to see that N(i) ≈ N 0 for all i meaning that the portion beyond k N 0 is non-negligible, and D(k) ≈ (k + 1)N 0 .

It is also possible to enhance our understanding of total deaths via an approximation involving the classical normal distribution. To this end, we work with the continuoustime counterparts: N c (t) for N(k) and D c (t) for D(k). That is, beginning with the infinitesimals over time interval dt

we integrate, substitute for d(t) and carry out a lengthy but straightforward calculation to obtain

is a scaling constant and ϕ(x) is the density function for a standard normal distribution; µ 0 . = θ 0 and σ 2 0 . = θ 0 d 0 / log(2), respectively, play the roles of mean and variance for a notional "time-to-peak" random variable. Thus, our approximation for total deaths is given by

where Φ(x) above is the cumulative distribution function for the standard normal distribution.

To provide an indication of the quality of the above approximation, in Figure 4 , we display D(k) along with D c (k) for six (θ 0 , d 0 ) combinations. It is clear that they are relatively close for the cases considered and we have found this to be true for a wide range of practical parameter values. Specifically, the maximum relative difference is approximately 11% for θ 0 , d 0 > 2 and 5% for θ 0 , d 0 > 7, but can be large if θ 0 and d 0 are both very small, which is not likely in practice. 

We now gain insight into the asymptotic total number of deaths by letting t tend to infinity and obtain

Of particular interest is the dependence on d 0 which could not be fully characterized for D(∞) previously. Thus, differentiating, we find that

where z 0 = µ 0 /σ 0 and A(z 0 ) = ϕ(−z 0 )/(1 − Φ(−z 0 )). Clearly, for θ 0 , d 0 > 0, the above derivative is negative if z 2 0 + A(z 0 )z 0 − 1 > 0. A straightforward numerical calculation shows that this holds true for z 0 > 0.84 which is equivalent to d 0 < 0.98θ 0 . Thus, D c (∞) decreases with d 0 until it reaches θ 0 , and, although this is an approximation to what happens for D(∞), it is nonetheless a useful insight into its behavior. In practice, we have found that θ 0 is larger than d 0 ; see the analysis of the two countries in Section 5 and note that this also true for a variety of other countries not shown. In summary, for practical purposes, one can view the total deaths as decreasing in d 0 .

Given a data set of actual deaths, N a (1), . . . , N a (n), for the purpose of estimating the parameters N 0 , d 0 , and θ 0 , using the resulting total death values D a (k), we minimize n k=1 (D a (k) − D(k)) 2 .

Since this objective function in nonlinear and non-convex in the parameters, we obtain good initial values by working with log-new-deaths. Then, a straightforward calculation leads to a classical linear least squares problem; i.e., minimizing

yieldsβ 0 ,β 1 , andβ 2 from which we obtain estimateŝ N 0 = 2β 0 ;d 0 = 1/(β 1 +β 2 );θ 0 = −(β 1 +β 2 )/(2β 2 ).

In some cases, these initial least squares estimates may be satisfactory solutions to the underlying problem in D a (k) and D(k). In other cases, further iterations are needed because they may be sensitive to noisy daily death data.

We now illustrate the application of our new model using historical data for year 2020 available in https://ourworldindata.org/coronavirus. To this end, consider death data for Brazil and Mexico beginning on March 28, 2020, a day on which Brazil had ninety-two total deaths and Mexico had only twelve, and ending on June 21, 2020. For both countries we estimated the parameter triple (N 0 ,d 0 ,θ 0 ) using the procedure described above. Then we compared model-based predictions with actual death numbers and made projections on the asymptotic value D(∞) obtained as k → ∞. The reader is reminded that an implicit assumption in our calculations is that the degree of mitigation u(k) is held constant; if a country prematurely relaxes measures such as social distancing, the model parameters should be recalibrated incorporating new data which comes to light.

In the case of Brazil, the estimated parameter values areN 0 = 28.0,d 0 = 6.68, and θ 0 = 70.4. The favorable performance of our model is depicted in Figure 5 where new deaths appear to be peaking around the time the model predicts. By summing up the daily totals, we readily obtain total deaths which increase sigmoidally to 49, 976 by the end of the observation period; interestingly, the asymptotic value D ∞ ≈ 74, 000 is about 50% higher than this. In the case of Mexico, the estimated parameter values areN 0 = 19.9, d 0 = 9.25, andθ 0 = 86.0. The fit to the data is shown in Figure 6 where, in contrast to Brazil, daily deaths had not quite peaked by the end of the observation period. Again by summing up new deaths, we obtain 20, 781 as the number of total deaths and note that the asymptotic value D(∞) ≈ 44, 000 is over 100% higher. 

This paper, part of the voluminous body of literature on epidemic modelling and control, differs from previous work using SIR-type models; i.e., we do not structure the dynamics based on many possible epidemiological considerations. Rather, we focus on a death doubling parameter d(k) and a peak day parameter θ(k), and rely on data to dynamically update the model as required. This approach to model identification is similar to those of [11] and [16] and, as previously mentioned, is motivated by the fact that each epidemic can present a vastly different array of challenges. For a new epidemic such as Covid-19, our view is that it may be premature to use highly structured model equations which rely on detailed epidemiological factors. As evidenced by our numerical experiments, our data-driven approach, with very few parameters to be estimated, appears to fit the data quite well; this is also true for simulations we conducted for many other countries not shown. Based on our work to date, two important directions immediately present themselves for future work: First, for the case when the mitigation level u(k) is no longer constant, it would be interest to study the evolution of deaths in an adaptive control context; i.e., as system parameters due to the arrival of new data, the controller u(k) is correspondingly adjusted. The control-theoretic setup in Figure 1 should rightfully be viewed in this more general context.

The second area for future research involves the formulation of an appropriate performance index. To provide some flavor as to the type of performance quantification issues which arise, for the current Covid-19 crisis, it is noted that societies around the world have been grappling with the following questions: How does the degree of mitigation get reflected in N(k) and D(k)? What tradeoffs is a society willing to accept between "lifestyle restrictions" and the level of deaths? Although existing literature includes "optimal control" formulations for epidemics, it is silent as to the detailed construction of the performance index. To illustrate what is meant by this, for the classical quadratic case

to say the least, it may be highly challenging to choose weights Q and R reflecting the tradeoffs which a society is willing to accept. In fact, it may prove to be the case that there are other cost functionals which are better suited for the study of pandemics.

A contribution to the mathematical theory of epidemics

Infectious Diseases of Humans: Dynamics and Control

The Mathematical Theory of Infectious Diseases and Its Applications

Networks, Crowds, and Markets: Reasoning About a Highly Connected World

A thousand and one epidemic models

Modelling the COVID-19 epidemic and implementation of populationwide interventions in Italy

A novel COVID-19 epidemiological model with explicit susceptible and asymptomatic isolation compartments reveals unexpected consequences of timing social distancing

Mathematical modelling of COVID-19 transmission and mitigation strategies in the population of Ontario, Canada

A mathematical model for simulating the phase-based transmissibility of a novel coronavirus

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)

Estimating initial epidemic growth rates

Transmissibility of 1918 pandemic influenza

Estimating the number of infections and the impact of nonpharmaceutical interventions on COVID-19 in 11 European countries

A modified SIR model for the Covid-19 contagion in Italy

Some epidemiological models with nonlinear incidence

Fitting dynamic models to epidemic outbreaks with quantified uncertainty: A primer for parameter uncertainty, identifiability, and forecasts

Optimum control of epidemics

Optimal vaccination schedules in a deterministic epidemic model

On the optimal control of a deterministic epidemic

Optimal control of deterministic epidemics

Stability analysis and optimal vaccination of an SIR epidemic model

Optimal control of vaccination and treatment for an SIR epidemiological model

An optimal control theory approach to nonpharmaceutical interventions

Optimal treatment of an SIR epidemic model with time delay

A new delay-SIR model for pulse vaccination

An optimal control problem for a spatiotemporal SIR model

A mathematical modeling with optimal control strategy of transmission of COVID-19 pandemic virus

On Fast Multi-Shot Epidemic Interventions for Post Lock-Down Mitigation: Implications for Simple Covid-19 Models

Optimal COVID-19 epidemic control until vaccine deployment