The COVID-19 situation and its immense toll on human lives has highlighted the enormous public health challenges associated with managing a pandemic. In the absence of a vaccine, there are primarily three control levers available for public health officials, namely, (a) contact restrictions, (b) testing, tracing and isolation, and (c) provisioning for additional medical capacity. Of these, contact restrictions via lockdowns and social distancing have emerged as the most powerful policy instrument especially in low and middle income countries that cannot afford to scale up testing or medical capacity. Choosing the optimal level of restrictions, however, has been highly non-trivial not only because it involves a complex trade-off between the yet to be understood COVID-19 disease impact and other socio-economic disruptions, but also because of the rapidly evolving situation on the ground.
Public health interventions related to the COVID-19 pandemic have largely been driven by scenariobased epidemiological forecasting studies [20, 17] . Current epidemiological models [40, 11, 7] incorporate spatio-temporal variations and predictive signals such as mobility to provide high fidelity forecasts. However, the decision making on contact restrictions is still fairly sub-optimal since it is based on the comparison of a few enumerated scenarios for a limited time horizon. Furthermore, forecasts based on a constant transmission rate (i.e., avg. #new people directly infected by an infectious person per time unit) convey the impression that the epidemic progression corresponds to a bell curve [40] regardless of empirical evidence to the contrary (see Fig 1) . Flattening the curve till herd immunity is seen as the only choice. Epidemiological analysis is often centred around determining the height and timing of the caseload peak as well as the time to attain herd immunity. Though highly valuable, this scenario-based decision-making approach leans towards a limited reactive role for public health agencies. In contrast, much less attention has been devoted to developing a mathematical control framework to support proactive decision-making based on the target disease & economic outcomes, and the state of the epidemic. Multiple studies [16, 29, 13, 41] point to the benefits of periodic lockdowns and staggered mobility among population groups, but these dynamic interventions are based on forecast simulations of a limited number of scenarios and are not adaptive in nature.
In this paper, we explore the problem of optimal adaptive control of transmission rate for a desired bound on infectious population. We focus on epidemiological models based on compartmental (SIR and SEIR) dynamics [28] because of their wide applicability, parsimonious & interpretable encoding of the disease dynamics, and amenability for data-driven calibration to yield accurate forecasts.
• We demonstrate that the SIR dynamics map to the well-known Lotka-Volterra (LV) system [8] on interpreting infectious patients as predators and susceptible contacts (i.e., the product of transmission rate and susceptible population) as the prey under specific conditions on the transmission rate. The resulting system (LVSIR) has a well-defined stable equilibrium point and an "energy" conservation property. It exhibits a bounded cyclic trend for active infections and a steady decline of the susceptible population. • We derive optimal control policy for transmission rate (CoSIR) using control-Lyapunov functions [45] based on the energy of the system, that is guaranteed to converge to the desired equilibrium, i.e., target infectious levels from any valid initial state. We also discuss extensions to compartmental model variants that involve an incubation period (e.g., delayed SIR, SEIR) as well as control of the infectious period that is influenced by testing and quarantine policy. • We propose a practical approximate implementation of the transmission rate control via discrete, but time-varying restriction levels. Simulation results demonstrate the efficacy of the approximate control in stabilizing infections and adaptability to perturbations from super-spreader events 2 .
The rest of the paper is organized as follows. Section 2 presents a formulation of the restriction control problem. Section 3 provides background on compartmental models, LV systems and relevant aspects of control theory. Sections 4, 5, 6, 7 present the SIR to LV system mapping, the transmission rate control mechanism, practical restriction control policy, and extensions respectively. Section 8 presents the concluding remarks.
Notation: x t and x(t) interchangeably denote the value of a variable x at time t. Time derivative of
t1 denotes the series of x as t varies from t 1 to t 2 .
During an epidemic, a key concern for public health officials is to determine the right level and schedule of contact restrictions that balances the disease and socioeconomic burdens. Strict lockdown conditions for a short time period suppress the infection levels, but infections tend to flare up again on easing restrictions unless the epidemic is completely wiped out. On the other hand, prolonged restrictions with no intermittent easing hinder economic activity and impose heavy costs on vulnerable population groups. Furthermore, progressive reduction in the susceptible population offers a chance for relaxation of restrictions that needs to be exploited.
Modelling the multi-faceted impact of contact restrictions requires accounting for region-specific cultural and economic constructs as well as the available medical capacity, a highly complex task. For tractability, we assume that the public health goal is to limit active infections to a certain target level determined via an independent impact analysis [2] . The controls available to the public health authorities can be viewed as multiple knobs that can be set to different levels (e.g., public transport at 50% occupancy). However, the need to communicate the policy to the general public and ensure compliance entails a simpler strategy centred around a few discrete restriction levels [4, 1] , (e.g. Table  1 ) and a preset schedule for a future time horizon, which is often longer than the intervals at which the epidemic observations are collected. For example, the infection levels might be monitored at a daily frequency, but the restriction guidelines (e.g., Level 2 on weekdays and Level 1 on weekends) might be chosen for a monthly period. Restriction Control. For a given region, let N, S curr , I curr be the total, current susceptible, and infectious populations respectively. Let I target avg be the target average infectious level. Let A be the set of restriction levels for which the transmission rate is known or can be estimated 3 Example. On Oct 1, 2020, a hypothetical city has a population of 13M of which 0.2M are currently infectious and 0.2M are post-infectious with the rest still susceptible. Assuming five restriction levels with transmission rates of [0.1, 0.2, 0.3, 0.4, 0.5], the objective is to figure out a restriction policy for the upcoming month so that the infectious count averages around I target avg = 150, 000.
Since our primary focus is on an analytical control framework, we make simplifying assumptions on the observability, (i.e., accurate estimates of the infectious population is possible via a mix of serosurveys and diagnosis tests) and the infection dynamics (region isolation, homogeneous interactions, negligible incubation period, and constant infectious period). Section 7 describes extensions when some of these assumptions are relaxed.
Our work builds on three research areas: (a) compartmental epidemiological models, (b) Lotka-Volterra systems, and (c) optimal control of non-linear dynamical systems.
Infectious diseases are commonly modeled using compartmental models where the population (N ) is divided into various compartments corresponding to different disease stages and the intercompartment transitions are governed by the model dynamics.
The SIR model [28, 21] is the simplest and most widely used one. The model comprises three compartments: Susceptible (S), Infectious (I) and Removed (R -includes immune & post-infectious persons) with the dynamics in Fig 2(b) . Here β is the rate of disease transmission from infectious to susceptible individuals, which largely depends on the contact restriction policy. γ is the inverse of the average infectious period, which depends on the testing and quarantine policy but is largely invariant when testing volumes are low [10] . In this model, the effective reproduction number (avg. #direct infections from each infection) is R ef f = βS γN . Existing restriction control approaches [44] are often guided by the principle of ensuring R ef f 1.
Certain infectious diseases have a significant incubation period when the individuals are infected but are not spreading the disease (non-infectious). The SEIR model [23] includes an additional E (exposed) compartment to model the incubation phase.
Lotka-Volterra (LV) systems [8, 48, 25] model the population dynamics of predator-prey interactions in a biological ecosystem. These models form a special case of Kolmogorov systems [22] that capture the evolution of a stochastic process over time. In a simple 2-species LV system, the population of prey (p) interacts with that of predator (q). The growth rate of prey depends on its reproduction rate (r) and the rate of consumption by predator (e). The change in predator population depends on the nourishment-based birth rate b and its death rate d. The system has two fixed points: (a) a saddle point that maps to extinction (p, q) = (0, 0), and (b) a stable equilibrium at (p, q) = ( d b , r e ). Typically, the system exhibits oscillations resulting in a closed phase plot that corresponds to conservation of an "energy" function. Fig 2(a) presents the dynamics of an LV system and the oscillations of the prey and predator populations. Due to the criticality of ecological population control, there has been considerable research on multiple variants of LV systems [22, 47] and their Hamiltonian dynamics [35] .
Optimal control of dynamical systems has rich connections to multiple fields [38] that deal with optimizing sequential decisions to maximize a desired objective such as reinforcement learning [15] , multi-armed bandits [12] , and stochastic control. Given a set of control variables, the optimal control policy describes the time derivatives of these variables that minimize the cost function and can be derived using Pontryagin's maximum principle [38] or the Hamilton-Jacobi-Bellman equations [37] .
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted November 13, 2020. ; Though there exist comprehensive techniques for control of linear dynamical systems, the control strategies for non-linear dynamical systems rely heavily on the existence of control-Lyapunov functions, which are typically identified using conservation laws of the associated physical systems. Once a suitable Lyapunov function is identified, there exist multiple control design strategies such as feedback linearization, backstepping and sliding mode control that are guaranteed to converge using Artstein's theorem [43] . In the case of the SIR model, a suitable Lyapunov function is not readily evident. On the other hand, Lyapunov stability and practical control strategies of LV systems have been extensively studied [34, 36, 32, 19 ].
Our primary goal is to solve the contact restriction control problem in Section 2. We choose to mainly focus on SIR dynamics because it captures the core disease spread mechanism that is an essential element of most epidemiological models. Since existing work on stability analysis of SIR models [5] does not address controllability, we first establish a connection between the SIR model and LV system, which is more amenable to principled control. In Section 5, we leverage the properties of the LV system to propose strategies for restriction control.
The problem of stabilizing infection levels assuming SIR dynamics has a direct analogy with population control in LV predator-prey systems where it is desirable to maintain the predator and prey population at certain target levels suitable for the ecosystem. Comparing the SIR and LV dynamics in Fig 2, we observe that the behaviour of the infectious people (I) is similar to that of the "predators" (q). There is an inflow (birth) βSI/N that depends on β as well as the current infectious & susceptible population and an outflow (death) γI from the I to the R compartment. However, the counterpart for the "prey" is not readily apparent.
An intuitive choice for "prey" is the "susceptible contacts" (i.e., the product of susceptible people and β, the #contacts of an infectious person per day) since this acts as "nourishment" to the predators and contributes to the inflow into the I compartment. Denoting the susceptible contacts by J = βS, we note that equivalence with the LV system requireṡ
(1) Hence, the transmission rate β has to folloẇ
This modified version of SIR model (LVSIR) maps to a Lotka-Volterra system. Comparing the model parameters, we note that the inverse of the infectious period (γ) corresponds to the predator death rate (d) and the inverse of population (1/N ) to the predator birth rate (b). We now analyze the behaviour of the LVSIR system.
Theorem 1 4 For the LVSIR model in Fig. 2(c) , the following holds true: 4 Please see Appendix A for the proofs.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 
where (x min , x max ) are defined as above, G(z) = γ(e z − z − 1) − w 0 , and F 1 (s), F 2 (s) are restrictions of F (s) = s + r log(1 − s r ) to +ve and -ve ranges. In general, T period has the form g(r, γ, w0 r ) with approximation via linearization yielding T period Table 2 . Similar to an LV system, the "energy" which corresponds to a weighted Itakura-Saito distance [26] between 6 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this this version posted November 13, 2020. ;
(I, J) and the equilibrium (I * , J * ) is conserved. The infectious population I (and susceptible contacts J) oscillates between [y min I * , y max I * ] (and [x min J * , x max J * ] ) during the entire period with an average value of I * (and J * ) while the susceptible population reduces steadily in a staircase-like fashion. The transmission rate β also exhibits periodic oscillations but the average steadily goes up to compensate for the reduction in the susceptible population. Fig 3(c-b) shows the phase plot and the variation of the key quantities during a single period with the four extreme points (South, East, North, West) marked explicitly. Fig 3(d) shows how T period [42] depends on r for different choices of w 0 /r. This mapping can be used to identify suitable values of r for practical policy making. For the special case where the initial state is at equilibrium, the system behaviour is steady as in Fig 3(a) .
We now consider the problem of controlling the transmission rate β for the LVSIR model (Fig 2(c) ) to nudge the infectious levels to a desired equilibrium. As discussed in Section 3, control of non-linear dynamical systems is typically achieved via Control-Lyapunov functions (CLFs). Our approach is to exploit the mapping from the SIR model to LV system and use CLFs derived from the "Lotka-Volterra energy" function w(J, I).
Given a dynamical systemż = f (z, u) with state vector z ∈ D ⊂ R n , control u ∈ R m , and equilibrium state z * = 0, a control-Lyapunov function (CLF) is a function V :
The CLF V (·) can be viewed as a generalized energy function withV (·) being a dissipation function. Artstein [6] proved that as long as there is a CLF, there exists a control u to ensure the reduction of energy at every non-equilibrium state and eventual convergence to the zero energy equilibrium.
,ż = f (z, u) = f 0 (z) + m i=1 f i (z)u i with state z ∈ D ⊂ R n , control u ∈ R m ,
Once a CLF is identified, it is relatively straightforward to design an appropriate control function u as described in [6, 46] . For our scenario, we rely on the conservation law of the LV system as well as the existing literature on its Lyapunov functions [39] .
Let z = ( J J * − 1, I I * − 1) so that the equilibrium z * = (0, 0). Let L(a 1 , a 2 ) : R 2 + → R be a continuously differentiable divergence such that | dL da1 | > 0 and L(a 1 , a 2 ) > 0, ∀a 1 = a 2 and L(a 1 , a 2 ) = 0 ⇐⇒ a 1 = a 2 . Then, the function V (z) = L(w, w target ) where w = w(J, I) can be used as a CLF. We focus on the case where w target = w * = 0. and propose a controlled SIR model (CoSIR).
Theorem 4 (CoSIR) For the SIR model, a proportional additive control on β, the convergence rate dependent on η. The design of the control also makes it robust to perturbations in the infectious population as the system re-calibrates β as appropriate.
The β-control policy (Eqn. 3) can be interpreted as follows. The first term β 2 I N corresponds to the relaxation possible due to the decreasing susceptible population while the second term (r − eI)β leads to oscillatory behavior, and the last term uβ = −η(t)β dL dw (J/J * − 1) ensures dissipation of energy and convergence to the equilibrium.
We now describe a practical solution to the public health restriction control problem in Section 2. Algorithm 1 outlines a holistic approach to obtain a restriction schedule using the optimal β-control in Theorem 4. There are four key steps.
Input Collection. Infection level targets (I target avg , I target max ), periodicity of the restriction schedule (T period ), decision horizon (T ) need to be determined based on a careful assessment of public health and socioeconomic considerations. Historical case counts and restrictions ([S t , I t , R t , a t ], [t] curr 0 ) also need to be collected to enable accurate optimization.
Data-driven Calibration. The next step is to use SIR-calibration methods [40, 11, 7] along with historical data to estimate a static γ, time-varying β, and the state of the epidemic (S curr , I curr , R curr ). The restriction level to transmission map, π, can be initially chosen from public health guidelines [1] and refined using the observed β for past restrictions in the region of interest.
Choosing COSIR Parameters. The free parameters of the COSIR model need to be chosen based on the control requirements. Algorithm 1 lists the updates derived from Theorem 2 with flexibility on the choice ofβ curr and η. Choosing the immediate transmission rate to beβ curr = J * Scurr = γN Scurr (equivalent to forcing effective reproduction number R ef f = 1) ensures a maximal reduction in the system "energy" and faster convergence to the desired equilibrium, but dampens fluctuations. However, fluctuations might be necessary for economic activity. When nearly steady infection levels are desired, I target avg I target max , then r = (2π) 2 γT 2 period and high η are appropriate.
Computing Optimal Restrictions. Determining the optimal policy can be split into two phases. The first involves estimating the ideal β control from Eqn. 4 while the second involves identifying the "closest" restriction level for the ideal β at each time step with "closeness" based on a suitable divergence such as the squared loss.
Employing Algorithm 1 for the scenario in Example 1 with a weekly periodicity, we obtain suitable parameters (Table 2 ) for the COSIR model. Fig 4 shows the behaviour with three variants of restriction policies (arbitrary, COSIR, and COSIR-APPROXIMATION based on the specified levels). In case of the arbitrary policy, the infections peak in narrow time intervals and become unmanageable while the ideal COSIR and even approximate control variants are able to achieve a steady rate of infections. This is true even when the infection levels are subject to sudden upward (t = 50) or downward perturbations (t = 100) in case of super-spreader events or sudden quarantine restrictions respectively. The β-control mechanism adapts and continues to push towards the equilibrium.
8 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. 
We now briefly describe some key extensions. Delayed SIR & SEIR Models. The SEIR model allows explicit modelling of the incubation period and is known to closely mimic the behaviour of the delayed SIR model [27] . When β follows Eqn. 2, the delayed SIR model readily maps to a delayed LV system with a non-preying growth period for the predators, which is a special case of the well-studied Wangersky-Cunningham systems [33, 47] . It can be shown that the modified delayed SIR system with a delay τ has the same equilibrium (J * , I * ) = (γN, r/e), exhibits (unbounded) oscillations and permits control of the form Eqn 3, where
There is a need for special handling when J approaches J * with the behavior depending on τ .
Testing & Isolation Policy. Testing, tracing and isolation also play a critical role in regulating the epidemic. In terms of SIR and SEIR dynamics, the net effect of aggressive testing is minimizing the infectious period or increasing γ [9, 30] . This is analogous to the culling of predators (infectious population) by increasing the death rate for which there already exist multiple control mechanisms [34] . In particular, choosing V (z) = L(w(J, I), w * ) as the control-Lyapunov function of interest, we obtain the control, γ(t) = γ 0 + ζ(t) dL dw ( I I * − 1), with ζ(t) > 0.
Online Learning. The restriction control problem can also be posed as a non-linear contextual bandit [12] formulation with the cumulative LV energy w(J, I) of the CoSIR model over a future horizon interpreted as the (negative) "reward". Here, the discrete restriction levels can be viewed 9 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. as the multiple arms of a bandit, the context includes the state of the epidemic, and the "reward" distribution is computed using the context and the observed transmission rate for the arms.
Our current work proposes an analytical framework for epidemic control with the intent of supporting an active goal-oriented public health response. The proposed framework relies on a mapping between SIR dynamics to Lotka-Volterra systems under specific conditions on the transmission rate. Given the vast literature on control of LV systems, this mapping can be leveraged to design new epidemic control techniques as well as extend current results to richer heterogeneous compartmental models and additional control variables (e.g., testing levels). Effective practical implementation of control requires further exploration of online and reinforcement learning variants [38] , addressing the limitations of SIR dynamics, and incorporation of additional signals such as mobility [3] . This effort also points to the feasibility of exploring control strategies for other macro social and physical systems, e.g., good and bad actor populations on social media.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted November 13, 2020. ; https://doi.org/10.1101/2020.11.10.20211995 doi: medRxiv preprint Denoting w b = (r + γ)D max , from Theorem 2(a), we note that w(J(t), I(t)) = w(J 0 ,
Given the nature of f (s), f (s) < c ⇒ s min < s < s max where (s min , s max ) are the finitevalued roots of f (s) = c. Hence x(t) and y(t) are both bounded on either side. Consequently, (x(t) − 1, y(t) − 1) is confined to a bounded rectangle and thus,
where can be directly expressed in terms of δ and vice versa. Hence, from Definition 1, z * = (1, 1) (or equivalently (J * , I * )) is a stable equilibrium.
When initial state (J 0 , I 0 ) is at equilibrium (J * , I * ) = (γN, r/e), we have (J(t), I(t)) = (γN, r/e), ∀t. Hence,Ṡ
Similarly,Ṙ = γI * ⇒ R(t) = R 0 + γI * t, 
Let z = log(x). Then,ẋ = e zż andẍ = e z (z +ż 2 ).
From 7, we have e z (z +ż 2 ) − e 2z (ż) 2 e z −γ(e z − 1)(e zż − re z ) = 0 ⇒z − γ(e z − 1)(ż − r) = 0.
Choosing s =ż ⇒ s =ẋ x = −r(y − 1).
Let w 0 = w(J 0 , I 0 ). Then, the trajectory corresponds to γ(x − log x − 1) + r(y − log y − 1) = w 0 ⇒ γ(e z − z − 1) + r(y − log y − 1) = w 0
where G(z) = γ(e z − z − 1) − w 0 and F (s) = s + r log(1 − s r ). Let F 1 (s), F 2 (s) be the restrictions of F (s) for the lower and upper parts of the phase plot. Then the time period for the lower section is given by dz.
. CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted November 13, 2020. ; https://doi.org/10.1101/2020.11.10.20211995 doi: medRxiv preprint When w 0 0, linearization is possible. Simplifying the trajectory F (s) = G(z) using the approximations e a = 1 + a + a 2 2 and log(1 − a) = a − a 2 2 , we have γ(e z − z − 1) − w 0 = s + r log(1 − s r )
Essentially, we have an elliptical curve with x, y following sinusoidal behavior with a period 2π √ rγ .
Proof of Theorem 2(d).
Assuming a continuous form for J, we observe thaṫ Idt (since I is periodic) = γI * T period (from above) 17 . CC-BY-NC-ND 4.0 International license It is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) preprint
The copyright holder for this this version posted November 13, 2020. ;