A Data-Driven Control-Theoretic Paradigm for Pandemic Mitigation with Application to Covid-19
aa r X i v : . [ q - b i o . P E ] A ug A Data-Driven Control-Theoretic Paradigmfor Pandemic Mitigation with Application toCovid-19
Kevin Burke ∗ B. Ross Barmish † August 17, 2020
Abstract
In this paper, we introduce a new control-theoretic paradigm for mitigating thespread of a virus. To this end, our discrete-time controller, aims to reduce the num-ber of new daily deaths, and consequently, the cumulative number of deaths. Incontrast to much of the existing literature, we do not rely on a potentially complexvirus transmission model whose equations must be customized to the “particulars”of the pandemic at hand. For new viruses such as Covid-19, the epidemiology driv-ing the modelling process may not be well known and model estimation with limiteddata may be unreliable. With this motivation in mind, the new paradigm describedhere is data-driven and, to a large extent, we avoid modelling difficulties by concen-trating on just two key quantities which are common to pandemics: the doublingtime , denoted by d ( k ) and the peak day denoted by θ ( k ). Our numerical studies todate suggest that our appealingly simple model can provide a reasonable fit to realdata. Given that time is of the essence during the ongoing global health crisis, theintent of this paper is to introduce this new paradigm to control practitioners anddescribe a number of new research directions suggested by our current results. Keywords.
Death data; Discrete time; Disease modelling; Doubling; Peaking
MSC 2010 subject classification: ∗ Corresponding author. University of Limerick, Ireland; [email protected] † Boston University and University of Wisconsin; [email protected]. Introduction
The main objective in this paper is to introduce a new control-theoretic paradigm formitigation of the spread of pandemic disease. In the sequel, the discrete-time controlvariable u ( k ) corresponds to the “degree of mitigation” which could be clinical or non-clinical in nature. The motivation for the paradigm to follow is simple to explain: Whena new pandemic emerges on the world scene, standard virus transmission models aretypically adopted in an effort to predict and control the uncertain future. However, withlimited epidemiological information available for a novel virus, existing models may beinappropriate; i.e., such a virus may have many features in its transmission and mitigation-response dynamics which are not well understood and not captured using virus modelsfrom the past. Accordingly, with this motivation in mind, data-driven modelling is thefocal point in this paper.In the existing literature, the majority of papers on virus transmission are based on theso-called “SIR” compartmental model[1] which, at a given point in time, places individualsinto 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 considerationssuch as death due to the disease, the latency time for infected individuals to becomeinfectious, 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 modelto the specifics of the situation at hand is non-trivial. This “customization” issue isepitomized 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 epidemiologicalconsiderations. In the early days of a pandemic, with limited data, not only is it difficultto build the structural model, but also to estimate its parameters. A particularly thornymodelling issue is that many infected individuals go undetected [10], in large part dueto asymptomatic or mild cases. This runs counter to many SIR modelling efforts whichassume that the detected cases are equivalent to the true cases when, in fact, the formerwill 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 startingpoint is to work with data in the form of death numbers, which are typically higher fidelityas these are much less likely to be undetected or under-reported; e.g., see [11]. Althoughat 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 inthe literature. In another recent paper [14], this issue is partially remedied by usingboth death and case data with some additional parameters to account for the differencebetween detected and true cases. In practice, the high variability in the detection processseems consistent with the wide uncertainty bands produced by [13] for the proportion ofinfected individuals in the population.In summary, while research efforts aimed at understanding virus transmission mech-anisms are valuable and important, the challenges in modelling will not be lost on the2eader. Therefore, as mentioned above, we focus entirely on the death process data andpropose a parsimonious three-parameter model. We expect that this approach will appealto the control community wherein epidemiology is not a core competency. That beingsaid, although we avoid transmission dynamics, we make use of two quantities which arecommon to epidemics over limited time frames: the so-called “doubling time” and “peakday” 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 andempirically-based phenomenological models in [11] and [16].In addition to the considerations above, in this paper, we also include the dynamics ofa controller reflecting the effect of mitigation measures aimed at reducing new deaths tozero. To the best of our knowledge, the earliest work bringing control-theoretic method-ologies to epidemiological modelling dates back to the 1970s; see [17, 18, 19]. Buildingon this, in [20] the theoretical properties of various mitigation techniques are considered,while others focussed more specifically on vaccination [21], the combination of vaccinationand treatment [22], non-clinical interventions such as social distancing, quarantining, andeducation [23], system delay [24, 25] and spatiotemporal effects [26]; more recently, theCovid-19 outbreak has been emphasized in [27, 28, 29]. Interestingly, all of these con-trol strategies assume underlying SIR-type models which, as previously discussed, presentvarious challenges in practice. In contrast to these approaches, a key feature of our data-driven approach is the error dynamics involving the comparison of model-predicted deathsto actual deaths. This enables both the model and the mitigation controller to be adaptedover time.The remainder of the paper is organised as follows: In Section 2 we describe the pre-liminaries involving the use of data and error dynamics at a high level. Subsequently, inSections 3 and 4, the details of our new model and its qualitative properties are providedfor the case of constant mitigation. In Section 5, we describe the process of model pa-rameter estimation and illustrate its use via numerical examples involving Covid-19 datafrom Brazil and Mexico, two of the countries among those with the highest of death ratesas of mid-2020. Finally, in Section 6, conclusions and directions for future research aredescribed 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 thecorresponding level of mitigation provided by the controller. In this first paper aimed atintroducing 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 asgovernment mandates on social distancing and the use of masks and therapeutics whereassmaller values might correspond to “relaxation” of the rules.Starting at k = 0, at a general level, one begins with an equation for new deaths N ( k + 1) = f ( D ( k ) , N ( k ) , u ( k ))3here the cumulative total deaths are naturally constrained to be D ( k ) = D ( k − N ( k ).We propose an attractive form for f ( D, N, u ) in Section 3, based on the doubling timeand peak day quantities from epidemiology.
Our control-theoretic paradigm begins with the acquisition of daily death data which iscollected to obtain estimates of the parameters describing the function f above. Once themodel is fixed, consistent with the tenets of receding horizon control, one can make predic-tions 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 zerolimit may not be required if the number of deaths is an acceptably small fraction of thepopulation size.As a pandemic unfolds over time, the model estimation should be updated periodicallyas new data is obtained. That is, letting N a ( k ) denote the “actual” number of new deathson day k , as depicted in Figure 1, we use the error e ( k ) . = N a ( k ) − N ( k ) to drive the updateof 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-weightthe distant past or smooth the noise using cumulative errors. However, such considerationsare beyond the scope of this article.Figure 1: Adaptive Mitigation Control Our analysis begins with the so-called doubling time , a widely reported metric used inpractice 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 thisterminology for cumulative total deaths. k + 1 is givenby N ( k + 1) = 2 /d N ( k ) so that N ( k + d ) = 2 N ( k ) as expected. In practice, d = d ( k )is time-varying for reasons such as immunity building up in the population, changes tomitigation strategies, and medical developments. Therefore, with initial values d (0) = d and N (0) = N , viewed as “baseline” quantities at the point from which we model, weobtain new deaths as N ( k + 1) = 2 /d ( k ) N ( k ) . In the equation above d ( k ) can be quite general in its functional form. We now specialisealong the lines described in Section 1. That is, we structure d ( k ) so that it is consistentwith many standard epidemic models which exhibit “peaking” behaviour. Thus, with d >
0, the count on new deaths N ( k ) climbs until some peak day θ ( k ) . = θ ( u ( k ))which may be time-varying due to changes in the level of mitigation characterized by thecontroller 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 d ( k ) = θ ( k ) θ ( k ) − k d . which increases as k approaches θ ( k ) > k ∗ , we have k < k ∗ <θ ( k ) and k > k ∗ > θ ( k ), then a single peak occurs at k = k ∗ = θ ( k ∗ ). More generallyhowever, θ ( k ) may be referred to as the “anticipated” peak day since it can move furtheror 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 isbeing held fixed over some extended time period. To this end, we consider the controller u ( k ) = u 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 ofthe efficacy of control measures in place. Subsequently, if the updated model parameterspredict 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 θ . = θ ( u ) . and we obtain the equation for new deaths N ( k + 1) = 2 /d ( k ) N ( k ) = 2 θ − kθ d N ( k )5hich is readily solved for k ≥
0. Indeed, expressing N ( k ) as a product followed bysummation of exponents yields N ( k ) = k − Y i =0 θ − kθ d N (0) = 2 (2 θ k − k θ d N . The appealingly simple solution for N ( k ) above lends itself to various insights about theevolution and qualitative behavior of deaths. First, it is evident by inspection that N ( k )is bell-shaped due to the negative k exponent, and this is clear from Figure 2. It is alsostraightforward to study the dependence of N ( k ) on θ and d by first noting that ∂ log N ( k ) ∂θ = k ( k −
1) log(2)2 θ d ,∂ log N ( k ) ∂d = k ( k − k N ) log(2)2 θ d where k N . = 2 θ + 1. For the most important case where d and θ are positive, we seethat N ( k ) increases in θ , decreases in d for k < k N , and increases in d for k > k N .It is apparent that k N represents the number of days until the pandemic reverts backto the early stage baseline level; i.e., N ( k N ) = N . However, in contrast to the earlystage, for k > k N , new deaths N ( k ) are now on a downward trajectory which in somesense can be viewed as signalling the end point of the pandemic. Finally, it is also ofinterest to characterize the number of deaths on the peak day, N ( θ ) = 2 θ d N which, clearly, increases in θ and decreases in d .The peak N ( θ ) above is important in that it correlates strongly with “anticipatedpressure” on hospitals; e.g., it is a predictor of stress on resources such as intensive careunits. In particular, given the concern that N ( θ ) exceeds some critical level N max , onecan easily use the equations above to study safety margins such as S ( d , θ ) = 100 (cid:18) − θ d N N max (cid:19) which indicates how far from N max , as a percentage, peak deaths will be. When S <
0, since this corresponds to N ( θ ) > N max , one might consider increasing the level ofmitigation by updating u ( k ) to reduce a potential crisis. In Figure 3, we display a realmof possible outcomes for this safety margin where, for example, with an initial doublingtime of 10 days, an anticipated peak of up to about 14 days, the ( d , θ ) point lies in thegreen-colored safety zone. 6
10 20 30 40 50 60
New Deaths
Day Number N u m b e r o f D ea t h s (14, 5)(14, 10) (28, 5)(28, 10) Figure 2: New Deaths Displayed for Four ( θ , d ) Cases with N = 1 −40−20020407 8 9 10 11 12 13 14810121416182022 Safety Margin Percentage Doubling Time P eak D ay Figure 3: Projected Safety Margin with N /N max = 0 . .2 Total Number of Deaths and Insights The total number of deaths, given by D ( k ) = k X i =0 N ( i ) , is not summable as a closed form solution. However, many properties of D ( k ), inheritedfrom N ( k ), are nevertheless clear: First, since each N ( i ) term is a point on a bell-shapedcurve, D ( k ) must be a nondecreasing sigmoidal function. Moreover, since N ( i ) is increas-ing in θ for all i , D ( k ) must also be increasing in θ . This makes the importance ofreducing θ by mitigation quite clear. That is, a reduced θ leads to a smaller peak whichhappens earlier, and, in turn, this lowers the total number of deaths. Furthermore, it isalso apparent that D ( k ) decreases with respect to d for k ≤ k N .For k > k N , the N ( i ) terms are increasing in d which makes the global behaviourof D ( k ) with respect to d non-trivial. However, recalling that k N can be viewed assignalling the end of the pandemic, for reasonable d , the terms entering D ( k ) beyond N ( k N ) will be quite small; i.e., the portion of D ( k ) which increases with respect to d will be small enough such that D ( k ) decreases with d . That being said, for very large d , it is easy to see that N ( i ) ≈ N for all i meaning that the portion beyond k N isnon-negligible, and D ( k ) ≈ ( k + 1) N . D ( k ) Using a Normal Distributions
It is also possible to enhance our understanding of total deaths via an approximationinvolving the classical normal distribution. To this end, we work with the continuous-time counterparts: N c ( t ) for N ( k ) and D c ( t ) for D ( k ). That is, beginning with theinfinitesimals over time interval dt dN c N c = log(2) d ( t ) , we integrate, substitute for d ( t ) and carry out a lengthy but straightforward calculationto obtain N c ( t ) = e R t d ( ζ ) dζ N = ϕ (cid:18) t − µ σ (cid:19) C N where C . = 2 θ d s πθ d log(2)is a scaling constant and ϕ ( x ) is the density function for a standard normal distribution; µ . = θ and σ . = θ d / log(2), respectively, play the roles of mean and variance for a8otional “time-to-peak” random variable. Thus, our approximation for total deaths isgiven by D c ( t ) = N + Z t N c ( ζ ) dζ = (cid:20) (cid:26) Φ (cid:18) t − µ σ (cid:19) − Φ (cid:18) − µ σ (cid:19)(cid:27) C (cid:21) N where Φ( x ) above is the cumulative distribution function for the standard normal distri-bution.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 ( θ , d ) combinations. It is clear that theyare relatively close for the cases considered and we have found this to be true for awide range of practical parameter values. Specifically, the maximum relative difference isapproximately 11% for θ , d > θ , d >
7, but can be large if θ and d areboth very small, which is not likely in practice. Total Deaths: Discrete Time (Solid); Continuous Time (Dash)
Day Number N u m b e r o f D ea t h s ( l og − sca l e ) (1, 5)(1, 10)(14, 5)(14, 10)(28, 5)(28, 10) Figure 4: Total Deaths Displayed for Six ( θ , d ) Cases with N = 1. D ( k ) We now gain insight into the asymptotic total number of deaths by letting t tend toinfinity and obtain D c ( ∞ ) = lim t →∞ D c ( t ) = (cid:20) (cid:26) − Φ (cid:18) − µ σ (cid:19)(cid:27) C (cid:21) N . Of particular interest is the dependence on d which could not be fully characterized for D ( ∞ ) previously. Thus, differentiating, we find that ∂ log D c ( ∞ ) ∂d = N − D c ( ∞ )2 d D c ( ∞ ) (cid:8) z + A ( z ) z − (cid:9) z = µ /σ and A ( z ) = ϕ ( − z ) / (1 − Φ( − z )). Clearly, for θ , d >
0, the abovederivative is negative if z + A ( z ) z − >
0. A straightforward numerical calculationshows that this holds true for z > .
84 which is equivalent to d < . θ . Thus, D c ( ∞ )decreases with d until it reaches θ , and, although this is an approximation to whathappens for D ( ∞ ), it is nonetheless a useful insight into its behavior. In practice, wehave found that θ is larger than d ; see the analysis of the two countries in Section 5and note that this also true for a variety of other countries not shown. In summary, forpractical purposes, one can view the total deaths as decreasing in d . Given a data set of actual deaths, N a (1) , . . . , N a ( n ), for the purpose of estimating theparameters N , d , and θ , using the resulting total death values D a ( k ), we minimize n X k =1 ( D a ( k ) − D ( k )) . Since this objective function in nonlinear and non-convex in the parameters, we obtaingood initial values by working with log-new-deaths. Then, a straightforward calculationleads to a classical linear least squares problem; i.e., minimizing n X k =1 (cid:18) log N a ( k )log 2 − β − β k − β k (cid:19) yields ˆ β , ˆ β , and ˆ β from which we obtain estimatesˆ N = 2 ˆ β ; ˆ d = 1 / ( ˆ β + ˆ β ); ˆ θ = − ( ˆ β + ˆ β ) / (2 ˆ β ) . In some cases, these initial least squares estimates may be satisfactory solutions to theunderlying problem in D a ( k ) and D ( k ). In other cases, further iterations are neededbecause they may be sensitive to noisy daily death data. We now illustrate the application of our new model using historical data for year 2020available in https://ourworldindata.org/coronavirus . To this end, consider deathdata for Brazil and Mexico beginning on March 28, 2020, a day on which Brazil hadninety-two total deaths and Mexico had only twelve, and ending on June 21, 2020. Forboth countries we estimated the parameter triple ( ˆ N , ˆ d , ˆ θ ) using the procedure describedabove. Then we compared model-based predictions with actual death numbers and madeprojections on the asymptotic value D ( ∞ ) obtained as k → ∞ . The reader is remindedthat an implicit assumption in our calculations is that the degree of mitigation u ( k ) isheld constant; if a country prematurely relaxes measures such as social distancing, themodel parameters should be recalibrated incorporating new data which comes to light.10n the case of Brazil, the estimated parameter values are ˆ N = 28 . , ˆ d = 6 .
68, andˆ θ = 70 .
4. The favorable performance of our model is depicted in Figure 5 where newdeaths appear to be peaking around the time the model predicts. By summing up thedaily totals, we readily obtain total deaths which increase sigmoidally to 49 ,
976 by the endof the observation period; interestingly, the asymptotic value D ∞ ≈ ,
000 is about 50%higher than this. In the case of Mexico, the estimated parameter values are ˆ N = 19 . , ˆ d = 9 .
25, and ˆ θ = 86 .
0. The fit to the data is shown in Figure 6 where, in contrast toBrazil, daily deaths had not quite peaked by the end of the observation period. Again bysumming up new deaths, we obtain 20 ,
781 as the number of total deaths and note thatthe asymptotic value D ( ∞ ) ≈ ,
000 is over 100% higher.
New Deaths (Green); 7-Day Rolling (Black); Prediction (Red)
Day Number Beginning March 28, 2020 N e w D ea t h s Figure 5: Daily Covid-19 Deaths in Brazil
New Deaths (Green); 7-Day Rolling (Black); Prediction (Red)
Day Number Beginning March 28, 2020 N e w D ea t h s Figure 6: Daily Covid-19 Deaths in Mexico11
Conclusion
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 dynamicsbased on many possible epidemiological considerations. Rather, we focus on a deathdoubling parameter d ( k ) and a peak day parameter θ ( k ), and rely on data to dynamicallyupdate the model as required. This approach to model identification is similar to thoseof [11] and [16] and, as previously mentioned, is motivated by the fact that each epidemiccan 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 whichrely on detailed epidemiological factors. As evidenced by our numerical experiments, ourdata-driven approach, with very few parameters to be estimated, appears to fit the dataquite well; this is also true for simulations we conducted for many other countries notshown.Based on our work to date, two important directions immediately present themselvesfor future work: First, for the case when the mitigation level u ( k ) is no longer constant, itwould be interest to study the evolution of deaths in an adaptive control context; i.e., assystem parameters due to the arrival of new data, the controller u ( k ) is correspondinglyadjusted. The control-theoretic setup in Figure 1 should rightfully be viewed in this moregeneral context.The second area for future research involves the formulation of an appropriate perfor-mance index. To provide some flavor as to the type of performance quantification issueswhich arise, for the current Covid-19 crisis, it is noted that societies around the worldhave 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 “op-timal control” formulations for epidemics, it is silent as to the detailed construction of theperformance index. To illustrate what is meant by this, for the classical quadratic case J ( u ) = Q N X k =1 N ( k ) + R N − X k =0 u ( k ) , to say the least, it may be highly challenging to choose weights Q and R reflecting thetradeoffs which a society is willing to accept. In fact, it may prove to be the case thatthere are other cost functionals which are better suited for the study of pandemics. References [1] W. O. Kermack and A. G. McKendrick, “A contribution to the mathematical theoryof epidemics,”
Proceedings of the Royal Society of London, Series A, vol. 115, pp. 700-721, 1927.[2] R. M. Anderson and R. M. May,
Infectious Diseases of Humans: Dynamics andControl,
Oxford University Press, 1991.123] N. T. J. Bailey,
The Mathematical Theory of Infectious Diseases and Its Applications ,Charles Griffin & Company Ltd, 1975.[4] D. Easley and J. Kleinberg,
Networks, Crowds, and Markets: Reasoning About aHighly Connected World,
Cambridge University Press, 2010.[5] H. Hethcote, “A thousand and one epidemic models,” in
Frontiers in mathematicalbiology, pp. 504-515, Springer, 1994.[6] G. Giordano, F. Blanchini, R. Bruno, P. Colaneri, A. Di Filippo, A. Di Matteo andM. Colaneri, “Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy,”
Nature Medicine, vol. 26, pp. 855-860, 2020.[7] J. Gevertz, J. Greene, C. H. S. Tapia and E. D. Sontag, “A novel COVID-19epidemiological model with explicit susceptible and asymptomatic isolation com-partments reveals unexpected consequences of timing social distancing,” medRxivpreprint, doi.org/10.1101/2020.05.11.20098335, 2020.[8] A.R. Tuite, D.N. Fisman and A.L. Greer, “Mathematical modelling of COVID-19 transmission and mitigation strategies in the population of Ontario, Canada,”
CMAJ, vol. 192, pp. E497-E505, 2020.[9] T.M. Chen, J. Rui, Q.P. Wang, Z.Y. Zhao, J.A. Cui and L. Yin, “A mathemat-ical model for simulating the phase-based transmissibility of a novel coronavirus,”
Infectious Diseases of Poverty, vol. 9, pp. 1-8, 2020.[10] R. Li, S. Pei, B. Chen, Y. Song, T. Zhang, W. Yang and J. Shaman, “Substan-tial undocumented infection facilitates the rapid dissemination of novel coronavirus(SARS-CoV-2),”
Science, vol. 368, pp. 489-493, 2020.[11] J. Ma, J. Dushoff, B. M. Bolker and D. J. Earn, “Estimating initial epidemic growthrates,”
Bulletin of Mathematical Biology, vol. 76, pp. 245-260, 2014.[12] C. E. Mills, J. M. Robins and M. Lipsitch, “Transmissibility of 1918 pandemic in-fluenza,”
Nature, vol. 432, pp. 904-206, 2004.[13] S. Flaxman, S. Mishra, A. Gandy et al., “Estimating the number of infections and theimpact of nonpharmaceutical interventions on COVID-19 in 11 European countries,”
Imperial College London COVID-19 Reports, doi.org/10.25561/77731.[14] G. C. Calafiore, C. Novara and C. Possieri, “A modified SIR model for the Covid-19contagion in Italy,” arXiv preprint,
Journal of Mathematical Biology, vol. 29, pp. 271-287, 1991.1316] G. Chowell, “Fitting dynamic models to epidemic outbreaks with quantified uncer-tainty: A primer for parameter uncertainty, identifiability, and forecasts,”
InfectiousDisease Modelling, vol. 2, pp. 379-398, 2017.[17] N. K. Gupta and R. E. Rink, “Optimum control of epidemics,”
Mathematical Bio-sciences, vol. 18, pp. 383-396, 1973.[18] H. W. Hethcote and P. Waltman, “Optimal vaccination schedules in a deterministicepidemic model,”
Mathematical Biosciences, vol. 18, pp. 365-381, 1973.[19] R. Morton and K. H. Wickwire, “On the optimal control of a deterministic epidemic,”
Advances in Applied Probability, vol. 6, pp. 622-635, 1974.[20] H. Behncke, “Optimal control of deterministic epidemics,”
Optimal Control Applica-tions and Methods, vol. 21, pp. 269-285, 2000.[21] G. Zaman, Y. H. Kang and I. H. Jung, “Stability analysis and optimal vaccinationof an SIR epidemic model,”
BioSystems, vol. 93, pp. 240-249, 2008.[22] T. T. Yusuf and F. Benyah, “Optimal control of vaccination and treatment for an SIRepidemiological model,”
World Journal of Modelling and Simulation, vol. 8, pp. 194-204, 2012.[23] F. Lin, K. Muthuraman, M. Lawley, “An optimal control theory approach to non-pharmaceutical interventions,” In
BMC infectious diseases, vol. 10, 2010.[24] G. Zaman, Y. H. Kang and I. H. Jung, “Optimal treatment of an SIR epidemic modelwith time delay,”
BioSystems, vol. 98, pp. 43-50, 2009.[25] C. Briat and E. I. Verriest, “A new delay-SIR model for pulse vaccination,” In
Biomedical Signal Processing and Control, vol. 4, pp. 272-277, 2009.[26] A. E. A. Laaroussi, M. Rachik and M. Elhia, “An optimal control problemfor a spatiotemporal SIR model,”
International Journal of Dynamics and Con-trol, vol. 6, pp. 384-397, 2018.[27] A. Kouidere, B. Khajji, A. El Bhih, O. Balatif and M. Rachik, “A mathematical mod-eling with optimal control strategy of transmission of COVID-19 pandemic virus,”
Communications in Mathematical Biologyanc Neuroscience, vol. 2020, ID. 24, 2020.[28] M. Bin, P. Cheung, E. Crisostomi, P. Ferraro, C. Myant, T. Parisini and R. Shorten,“On Fast Multi-Shot Epidemic Interventions for Post Lock-Down Mitigation: Impli-cations for Simple Covid-19 Models,” arXiv preprint, medRxivpreprint,medRxivpreprint,