A Closed-Loop Framework for Inference, Prediction and Control of SIR Epidemics on Networks
AA Closed-Loop Framework for Inference, Prediction and Control ofSIR Epidemics on Networks
Ashish R. Hota, Jaydeep Godbole, Pradhuman Bhariya and Philip E. Par´e ∗ Abstract
Motivated by the ongoing pandemic COVID-19, we propose a closed-loop framework thatcombines inference from testing data, learning the parameters of the dynamics and optimalresource allocation for controlling the spread of the susceptible-infected-recovered (SIR) epidemicon networks. Our framework incorporates several key factors present in testing data, such ashigh risk individuals are more likely to undergo testing and infected individuals potentially actas asymptomatic carriers of the disease. We then present two tractable optimization problems toevaluate the trade-off between controlling the growth-rate of the epidemic and the cost of non-pharmaceutical interventions (NPIs). Our results provide compelling insights for policy-makers,including the significance of early testing and the emergence of a second wave of infections ifNPIs are prematurely withdrawn.
Mathematical modeling of infectious diseases that spread through the human population has along history. A plethora of models that capture the dynamic evolution of epidemics have beenproposed (see [Hethcote, 2000, Pastor-Satorras et al., 2015] for detailed surveys). However, thereare two broad classes of models: one where a person who has recovered from the disease is immunefrom re-infection (at least for a few months or years) and second where there is a possibility ofre-infection.One of the most fundamental mathematical models pertaining to the first class of epidemics isthe susceptible-infected-recovered (SIR) epidemic where individuals can be in one of three possiblecompartments or states: susceptible, infected or recovered. Individuals who are susceptible getpotentially infected by coming into contact with infected neighbors, while infected individualsrecover at a certain rate. While several other models with additional compartments (such asdeath, quarantine, exposed, etc) have been proposed, the SIR epidemic remains one of the mostfundamental models. Early work on (SIR) epidemics focused on the homogeneous populationsetting, while more recent works model the interactions between individuals via an underlyingnetwork. Detailed overviews of different epidemic models over networks are given in [Barrat et al.,2008, Draief and Massouli, 2010, Pastor-Satorras et al., 2015, Nowzari et al., 2016, Mei et al., 2017].Since there is no strong evidence of re-infection in the ongoing COVID-19 pandemic, it may beassumed to fall under the first class of epidemics [Ota, 2020]. Consequently, we focus on estimationand control of SIR epidemics on networks in this work. Specifically, we • highlight the challenges associated with computing optimal NPIs to control the SIR epidemicon networks, ∗ Ashish R. Hota, Jaydeep Godbole and Pradhuman Bhariya are with the Department of Electrical Engineering,IIT Kharagpur, India. Philip E. Par´e is with the School of Electrical and Computer Engineering, Purdue University,USA. E-mail: [email protected], [email protected], [email protected], [email protected]. a r X i v : . [ phy s i c s . s o c - ph ] J un identify key characteristics inherent in testing data that have not been systematically capturedby the existing work on inference and prediction of epidemics, • rigorously establish the behavior of the discrete-time SIR epidemic on networks, • propose a stochastic non-linear observer model that relates testing data with the underlyingepidemic states, • propose a closed-loop framework that integrates inference, learning and optimal allocation ofNPIs, and • uncover compelling insights on the behavior of the SIR epidemic under different NPIs viaextensive empirical evaluations.Nevertheless, the proposed estimation and control approaches can naturally be extended to a largeclass of models with more compartments provided that there is no possibility of re-infection. In order to provide better clarity regarding the gaps in the existing literature and how our pro-posed framework fills those gaps, we first introduce the continuous-time SIR epidemic dynamics onnetworks [Mei et al., 2017]. Let G = ( V, E ) be a network or graph where V is the set of nodes with | V | = n and E ⊆ V × V is the set of edges. We consider a large-population regime where each noderepresents a sub-population (such as a city or a county/district or a state) as opposed to a singleindividual. We denote the size of sub-population i as N i .We denote by β ij ∈ R ≥ the rate at which the infection can spread through the edge ( v j , v i ) ∈ E ,and by γ i ∈ R > the rate at which an individual in sub-population i recovers from the infection.If two nodes v i and v j are not neighbors, then β ij = 0. We assume that β ii (cid:54) = 0 for everynode v i since the individuals inside sub-population i come in contact with each other. We define V i := { j ∈ V | β ij (cid:54) = 0 } to be the set of neighbors of node v i . In the literature, the notation a ij isoften used to denote the weight or contact pattern between nodes i and j and β i is used to denotethe rate at which node i is infected by its infected neighbors. The formulations are equivalent if wedefine β ij = β i a ij .The proportions of the sub-population at node v i that are susceptible, infected and recoveredat time t are denoted by s i ( t ) , x i ( t ) and r i ( t ), respectively. Accordingly, we have s i ( t ) , x i ( t ) , r i ( t ) ∈ [0 ,
1] and s i ( t ) + x i ( t ) + r i ( t ) = 1 for all t and v i ∈ V . The deterministic continuous-time evolutionof the SIR epidemic is given by˙ s i ( t ) = − s i ( t ) (cid:2) n (cid:88) j =1 β ij x j ( t ) (cid:3) , (1a)˙ x i ( t ) = s i ( t ) (cid:2) n (cid:88) j =1 β ij x j ( t ) (cid:3) − γ i x i ( t ) , (1b)˙ r i ( t ) = γ i x i ( t ) , ∀ i ∈ [ n ] := { , , . . . , n } . (1c)Note that ˙ s i ( t ) + ˙ x i ( t ) + ˙ r i ( t ) = 0. The evolution of the proportion or fraction of infected nodescan be written in a compact form as˙ x ( t ) = (cid:2) diag ( s ( t )) B − diag ( γ ) (cid:3) x ( t ) , (2)2here x ( t ) ∈ R n , B is the matrix with ( i, j )-th entry β ij and diag ( γ ) is a diagonal matrix withdiagonal entries being the entries of the vector γ . As analyzed in [Mei et al., 2017], the trajectoryof the infected population depends critically on the largest eigenvalue of the matrix (cid:2) diag ( s ( t )) B − diag ( γ ) (cid:3) , denoted by λ max ( t ). If λ max (0) >
0, then x ( t ) initially shows an exponential growth.Eventually, there exists a ¯ t > λ max (¯ t ) <
0, and the fraction infected monotonically decreasesfor t ≥ ¯ t . However, the proportion of the susceptible (recovered) population is monotonicallydecreasing (increasing) with t .Many infectious diseases, including COVID-19, show the above characteristic and such diseases,if unchecked, can potentially infect millions of people throughout the world in the span of severalweeks. In the absence of appropriate medicine and vaccines at the early stages of an outbreak (as isthe case for COVID-19) non-pharmaceutical intervention (NPI) strategies must be deployed. Theprimary NPIs to control the spread of such epidemics include • reducing the interaction between nodes (for example, by restricting travel and imposing lock-down measures), and • deploying more resources (in terms of healthcare personnel, dedicated hospitals and medicalequipment) thereby increasing the rate at which infected individuals get cured.The former corresponds to reducing the β ij parameters and the latter corresponds to increasing thecuring rates ( γ i ’s). Both of these interventions are costly; the former with a significant economiccost. Furthermore, different nodes have different degrees of epidemic outbreak, and require differentdegrees of interventions. Following early works in [Wan et al., 2007, 2008, Preciado et al., 2014], the existing literature onresource allocation for controlling the spread of epidemics has primarily focused the SIS epidemicon networks. In contrast with the SIR epidemic, under the SIS epidemic an individual after re-covery can potentially be infected again if she/he comes in contact with other infected individuals.Nevertheless, for the SIS epidemic, there is a simple spectral condition which characterizes whetherthe disease persists in the population or the fraction of infected nodes/populations decays quicklyto zero. It was shown in [Preciado et al., 2014] that the problem of optimal resource allocation toensure that the disease is eradicated is an instance of a geometric program (GP) [Boyd et al., 2007]which can be solved efficiently. This approach has been extended in several directions, such as toaccount for dynamic networks [Ogura and Preciado, 2017, Ogura et al., 2019], uncertainty in thenetwork structure [Han et al., 2015], distributed algorithms [Ram´ırez-Llanos and Mart´ınez, 2018,Mai et al., 2018], among others.
However, the above optimization problems are solved off-line anddo not use feedback to adapt the solution as the epidemic spreads. In contrast, investigations of optimal resource allocation to contain the spread of the SIR epi-demic have been few with [Ogura and Preciado, 2016] being an exception. It is shown in [Oguraand Preciado, 2016] that a GP can be formulated to minimize the expected cumulative number ofpeople who get infected when λ max (0) <
0. However, this condition is restrictive since for mostepidemics the fraction of infected individuals shows an initial exponential increase before declining.Optimal resource allocation for the SIR epidemic is particularly challenging because Online approaches based on optimal control theory have been studied in [Eshghi et al., 2014, Zaman et al., 2008]to compute vaccination strategies, which are different from NPIs considered in this work. there is no known tractable characterization of the eventual number of recovered individualsfor the SIR epidemic which can be minimized in an off-line manner, and • the instantaneous growth rate λ max ( t ) depends on the fraction of the susceptible populationat each node which is time-varying and may not be accurately known.In this paper, we investigate a potential approach to minimize λ max ( t ) in an online manner.As mentioned above, this requires knowledge of the proportion of the susceptible sub-populationat each node of the network. As a result, the current proportions of susceptible, infected, andrecovered sub-populations at different nodes and the parameters that govern the dynamics of theepidemic (such as infection and curing rates) need to be learned from testing data (number of testscarried out, number of confirmed cases and number of recoveries) that is made available every dayby different jurisdictions.Earlier work has primarily focused on learning the parameters of the epidemic dynamics. Inparticular, [Par´e et al., 2020] presents a data-driven framework for learning the parameters of net-worked SIS epidemic dynamics. However, literature on SIR epidemic models have mostly focusedon the scalar dynamics (without any network structure). In [Chen and Qiu, 2020], a least squaresbased parameter identification was carried out assuming that the proportion of infected and recov-ered sub-population and the daily change in the above (i.e., the state information) is proportionalto the respective fractions in the testing data. An analogous assumption was made in a recent work[Casella, 2020] on COVID-19 as well. In [Osthus et al., 2017, Song et al., 2020], Bayesian MarkovChain Monte Carlo (MCMC) techniques were used for estimating the states and parameters. Anoverview of the above approaches is discussed in [Abadie et al., 2020].However, the above models do not consider the networked SIR epidemic dynamics, and do notcapture the following characteristics inherent in testing data. • Real testing strategy is not uniform. Due to limited testing capacity, high-risk (symptomaticwith travel history) individuals are more likely to get tested [Cohen and Kupferschmidt, 2020]. • For some diseases, such as COVID-19, individuals often show symptoms a few days afterbecoming infected (while they continue to infect others despite being asymptomatic).
In light of the above research gaps, we propose a closed-loop framework that integrates inferringthe state of the epidemic from testing data, learning the parameters of the epidemic dynamicsand optimal resource allocation for the SIR epidemic on networks. The schematic of the proposedscheme is shown in Fig. 1.In Section 2, we first introduce a discrete-time counterpart of the SIR epidemic dynamics statedin (1), and provide assumptions required for the model to be well defined. We then derive severalproperties of the state trajectories under the discrete-time dynamics, including expanding the ideaof the reproduction number to the networked SIR model and showing at least linear convergenceto equilibria where no one is infected after some time.As discussed above, one of the key requirements towards controlling the evolution of the infectedpopulation is to infer the current fraction of susceptible population from testing data. We denotethe testing data at node i for a given time k as Ω i ( k ) := ( z i ( k ) , c i ( k ) , d i ( k )), where z i ( k ) denotes Estimating the final size of the recovered population is challenging, and existing approaches rely on approxima-tions that are not amenable for tractable optimization [Miller, 2012]. In epidemiology, the reproduction number is the number of infections one infection generates on average over thecourse of its infectious period. If less than one, the virus quickly dies out. ptimizationSolver Discrete-TimeSIR DynamicsParameterIdentification Data fromTestingInference B ∗ k , γ ∗ k θ ( k )Ω( k ) (cid:98) θ ( k − τ ) (cid:98) θ ( k ) Figure 1: Schematic of the proposed framework. Here θ ( k ) = ( s ( k ) , x ( k ) , r ( k )) denotes the epidemicstate at time k , Ω( k ) denotes the testing data, (cid:98) θ ( k ) denotes the inferred epidemic states and τ denotes the delay factor as explained in Section 1.3. the number of tests carried out, c i ( k ) denotes the number of confirmed cases and d i ( k ) denotes thenumber of recoveries. While prior work has assumed that c i ( k ) z i ( k ) is proportional to the (change in)fraction of infected population ( x i ( k )) [Chen and Qiu, 2020, Casella, 2020], we incorporate the factthat individuals at a higher risk of being infected are more likely to undergo testing than healthyindividuals. Furthermore, in order to capture the fact that for some diseases (such as COVID-19) an infected individual can remain asymptomatic for several days while being contagious, weassume that testing data is reflective of the epidemic state at a past time rather than the currenttime. We denote this lag or delay by τ . In practice, τ is different for each individual, and couldbe modeled as a random variable. However, as a first step, we assume τ to be deterministicand homogeneous across the population. In Section 3, we propose a stochastic nonlinear observermodel that relates the testing data with the underlying states by incorporating the two factorsdiscussed above. By leveraging the above model, we discuss how to infer the epidemic states (cid:98) θ ( k − τ ) = ( (cid:98) s i ( k − τ ) , (cid:98) x i ( k − τ ) , (cid:98) r i ( k − τ )) from testing data at time k .The delay factor τ necessitates estimating the current state of the epidemic, (cid:98) θ ( k ) :=( (cid:98) s ( k ) , (cid:98) x ( k ) , (cid:98) r ( k )), from the inferred values (cid:98) θ ( k − τ ). If the parameters of the epidemic dynam-ics ( β ij ( κ )’s and γ i ( κ )’s) are known for κ ∈ [ k − τ, k ], then the current state can be predictedfrom the dynamics. In the initial stages (before any interventions are deployed), these parametersmay be learnt from the available data. Eventually, once the optimal contact and curing rates aredeployed, they may be used for predicting the current state. We present a least-square formulationto identify the virus spread parameters in Section 4 in order to achieve this goal.In Section 5, we formulate two complementary GPs to evaluate the trade-off between minimizingthe cost of NPIs and the growth-rate (analogous expression of λ max ( t ) for the discrete-time model).The optimal solutions are then deployed, and at the next iteration, new testing results are usedto update the parameters and the predicted states, and the process repeats over time in a mannersimilar to that of receding horizon control. Finally, in Section 6, we illustrate the performance of the proposed inference, prediction andNPI allocation in controlling the spread of the SIR epidemic on networks. Our results highlightthe importance of early testing for accurate estimation and control performance, and the risk ofa second wave of infection if NPIs are prematurely withdrawn. We now present the discrete-timedynamics, inference, parameter identification and optimization formulations in the following foursections, respectively. Similar approaches have been investigated for the class of SIS epidemics in [Watkins et al., 2019]. However, tothe best of our knowledge, they have not been studied in the context of SIR epidemics. Discrete-Time SIR Epidemic Dynamics on Networks
The continuous-time dynamics in (1) is often viewed as a mean-field approximation of a continuous-time Markov chain model of the evolution of the disease. However, since the (testing) data is oftenavailable in discrete points of time, we develop the proposed inference, learning and control tasksfor a discrete-time version of the above dynamics motivated by a similar approach proposed in [Par´eet al., 2020] for SIS epidemics. We now state the discrete-time SIR epidemic dynamics obtained viaEuler discretization of (1). For a small enough sampling time h >
0, the deterministic discrete-timeevolution of the SIR epidemic is approximately given by s i ( k + 1) = s i ( k ) + h (cid:2) − s i ( k ) n (cid:88) j =1 β ij x j ( k ) (cid:3) , (3a) x i ( k + 1) = x i ( k ) + h (cid:2) s i ( k ) n (cid:88) j =1 β ij x j ( k ) − γ i x i ( k ) (cid:3) , (3b) r i ( k + 1) = r i ( k ) + hγ i x i ( k ) . (3c)The initial conditions need to be specified such that s i (0) , x i (0) , r i (0) ∈ [0 ,
1] and s i (0) + x i (0) + r i (0) = 1 for every node i . In vector form the model becomes s ( k + 1) = s ( k ) − h diag ( s ( k )) Bx ( k ) , (4a) x ( k + 1) = x ( k ) + h diag ( s ( k )) Bx ( k ) − h diag ( γ ) x ( k )= [ I n + h diag ( s ( k )) B − h diag ( γ )] x ( k ) =: A k x ( k ) , (4b) r ( k + 1) = n − s ( k + 1) − x ( k + 1) , (4c)where n is the vector of dimension n with all entries equal to 1. Note that we have defined A k := I n + h diag ( s k ) B − h diag ( γ ). We now assume the following on the parameters such that thedynamics is well behaved. Assumption 1.
For all i ∈ [ n ] , we have < hγ i ≤ and h (cid:80) nj =1 β ij < . Furthermore, the matrix B is irreducible. Note that Assumption 1 is satisfied when the sampling parameter h is chosen to be sufficientlysmall. If h is not sufficiently small (i.e., sampling is infrequent), then it is possible for the statesto become negative or exceed 1, both of which are incompatible with their physical interpretation.The assumption also implies that A k is an irreducible non-negative matrix. Therefore, by thePerron-Frobenius Theorem for irreducible non-negative matrices [Varga, 2000, Theorem 2.7 andLemma 2.4], A k has a positive real eigenvalue equal to its spectral radius, which, we denote by λ A k max .We have the following result on the behavior of the discrete-time dynamics under the aboveassumption. These are similar to the behavior of the continuous-time dynamics [Mei et al., 2017],but to the best of our knowledge, have not been formally proven in the literature. Our result alsostrengthens some of the observations in [Mei et al., 2017]. Theorem 1.
Consider the model in (4) under Assumption 1. Suppose s i (0) , x i (0) , r i (0) ∈ [0 , , s i (0) + x i (0) + r i (0) = 1 for all i ∈ [ n ] and x i (0) > for some i . Then, for all k ≥ and i ∈ [ n ] ,1) s i ( k ) , x i ( k ) , r i ( k ) ∈ [0 , and s i ( k ) + x i ( k ) + r i ( k ) = 1 ,2) s i ( k + 1) ≤ s i ( k ) , ) lim k →∞ x i ( k ) = 0 for i ∈ [ n ] ,4) λ A k max is monotonically decreasing as a function of k ,5) there exists ¯ k such that λ A k max < for all k ≥ ¯ k , and6) there exists ¯ k , such that x i ( k ) converges linearly to for all k ≥ ¯ k , i ∈ [ n ] . Further, if lim k →∞ s i ( k ) = 0 and hγ i = 1 for all i ∈ [ n ] , there exists ¯¯ k , such that x i ( k ) convergessuperlinearly to for all k ≥ ¯¯ k .Proof. See Appendix A.Theorem 1 shows that the model in (4) is well-defined, the susceptible proportions and thegrowth rate decrease monotonically over time, the growth rate will eventually be less than 1, andthe infected proportion of the population will go to 0, in at least linear time for large enough k . Remark 1.
The largest eigenvalue λ A k max is a generalization of the reproduction number to thenetworked epidemic setting, that is, if λ A k max < the virus quickly dies out. Our result rigorouslyproves the observation that for epidemics that follow an SIR-type dynamic, such as COVID-19, thereproduction number eventually falls below . We use these results as the framework for the control techniques presented in Section 5, whichalso require inferring the states from the testing data (Section 3) and estimating the spread param-eters for forecasting the states (Section 4).
As discussed earlier, part of the challenge in controlling the spread of infectious diseases such asCOVID-19 is that the prevalence (i.e., the underlying state) of the disease in a given populationis highly uncertain. Individuals need to be tested in order to determine if they are infected bythe disease or pathogen under consideration. At early stages of the epidemic, many countries andregions do not have enough capacity to test a large number of people since their testing kits arelimited. As a result, testing is conducted on individuals who show symptoms (such as fever andshortness of breath associated with COVID-19). However, drawing inference about the underlyingspread of the epidemic from such testing data is not straightforward since (i) a large fraction ofinfected and contagious individuals never show any symptoms, and (ii) similar symptoms are alsoexhibited by patients who suffer from other related illnesses [Hu et al., 2020]. Furthermore, testingdata on a given day reveals individuals who became infected a few days earlier (as opposed toinformation about the new infections on that day).In this section, we focus on inferring the epidemic state ( s ( k ) , x ( k ) , r ( k )) from testing data.Specifically, we develop a nonlinear observer model that relates the observations with the underlyingepidemic states following a Bayesian approach. We then leverage the proposed model for stateinference. We denote the inferred states with the symbol (cid:98) · .Recall that the testing data at time k at node v i is denoted by Ω i ( k ) and consists of the numberof tests carried out ( z i ( k )), the number of confirmed ( c i ( k )) and removed (including both recoveriesand deaths) ( d i ( k )) cases. We denote the cumulative number of confirmed and removed cases by C i ( k ) := (cid:80) kl =0 c i ( l ) and D i ( k ) := (cid:80) kl =0 d i ( l ), respectively. Accordingly, the number of known oractive cases is given by A i ( k ) := C i ( k ) − D i ( k ) and the daily change in the number of known/activecases is given by a i ( k ) = c i ( k ) − d i ( k ). 7 .1 Nonlinear Observer Model of Testing Data Note from the above discussion that testing enables us to observe the infection states of tested andconfirmed individuals. In particular, the known active cases are analogous to the proportion ofinfected individuals. As a result, it is reasonable to treat the proportion of new infections at time k (i.e., c i ( k ) z i ( k ) ) as representative of the fraction of new infections at node v i . However, as discussedearlier, the confirmed cases at time k are often found to have caught the infection several daysprior to being tested and with the delay denoted by τ ≥
1. Therefore, we assume c i ( k ) z i ( k ) to berepresentative of the decrease in the proportion of susceptible individuals at time k − τ , i.e., thequantity − ∆ s i ( k − τ ) = ∆ x i ( k − τ ) + ∆ r i ( k − τ ). In a departure from prior work which assumes c i ( k ) z i ( k ) to be proportional to x i ( k ) and/or − ∆ s i ( k ), we propose a Bayesian framework to model thefact that testing strategies are not uniform, and formally relate c i ( k ) z i ( k ) to ∆ s i ( k − τ ).We assume that a proportion h i ( k ) of the population at node v i belongs to a high risk categorywhile l i ( k ) := 1 − h i ( k ) comes under a low risk category. Individuals may be deemed high risk ifthey exhibit symptoms associated with the disease, are known contacts of confirmed cases, have atravel history in affected regions, or any combinations thereof.Now, let H i ( k ) be a random variable with H i ( k ) = 1 (resp. H i ( k ) = 0) if a randomly chosenindividual belongs to the high risk (resp. low risk) category. Accordingly, P ( H i ( k ) = 1) = h i ( k ).We also define a { , } -valued random variable DX i ( k ) with DX i ( k ) = 1 (resp. DX i ( k ) = 0) if arandomly chosen individual became infected at time k . We now introduce the following notationto denote certain conditional probabilities of interest. In particular, we define p hx,i ( k ) := P ( H i ( k ) = 1 | DX i ( k − τ ) = 1) = ⇒ P ( H i ( k ) = 0 | DX i ( k − τ ) = 1) = 1 − p hx,i ( k ) ,p hh,i ( k ) := P ( H i ( k ) = 1 | DX i ( k − τ ) = 0) = ⇒ P ( H i ( k ) = 0 | DX i ( k − τ ) = 1) = 1 − p hh,i ( k ) , where p hx,i ( k ) , p hh,i ( k ) ∈ [0 , p hx,i ( k ) captures the proportion of high risk in-dividuals (i.e., those who show symptoms) who became infected at k − τ and depends on thecharacteristics of the epidemic and the population (e.g., age, level of immunity against the disease,prevalence of comorbidity) at node v i . The parameter p hh,i ( k ) captures the proportion of healthyindividuals who belong to the high risk category (i.e., they show similar symptoms or have travelhistory), but did not become infected τ time steps earlier.With the above notation in place, we now apply Bayes’ law to compute the probability of arandomly chosen high risk individual being infected τ time steps earlier. Specifically, we compute P ( DX i ( k − τ ) = 1 | H i ( k ) = 1) = P ( H i ( k ) = 1 | DX i ( k − τ ) = 1) P ( DX i ( k − τ ) = 1) P ( H i ( k ) = 1)= p hx,i ( k )( − ∆ s i ( k − τ )) p hx,i ( k )( − ∆ s i ( k − τ )) + p hh,i ( k )(1 − ( − ∆ s i ( k − τ ))) =: p xh,i ( k ) . (5)Similarly, the probability of a randomly chosen low risk individual being infected τ time stepsearlier is P ( DX i ( k − τ ) = 1 | H i ( k ) = 0) = P ( H i ( k ) = 0 | DX i ( k − τ ) = 1) P ( DX i ( k − τ ) = 1) P ( H i ( k ) = 0)= (1 − p hx,i ( k ))( − ∆ s i ( k − τ ))(1 − p hx,i ( k ))( − ∆ s i ( k − τ )) + (1 − p hh,i ( k ))(1 − ( − ∆ s i ( k − τ ))) =: p xl,i ( k ) . (6)Recall that − ∆ s i ( k ) = s i ( k − − s i ( k ) = hs i ( k − (cid:80) nj =1 β ij x j ( k − ≥
0. Furthermore, both p xh,i ( k ) and p xl,i ( k ) are monotonically increasing in − ∆ s i ( k ).8e now relate the number of tests and confirmed positive cases with the underlying states. Attime k , the authorities at node v i decide to carry out z h,i ( k ) ∈ [0 , z i ( k )] number of tests on highrisk individuals and z i ( k ) − z h,i ( k ) number tests on low risk individuals. Assuming that the testingis accurate (with false positive and false negative rates being 0), the confirmed cases among thosetested is the sum of the number of confirmed cases among tested high risk individuals and amongtested low risk individuals. Accordingly, we model c i ( k ) ∼ Bin ( z h,i ( k ) , p xh,i ( k )) + Bin ( z i ( k ) − z h,i ( k ) , p xl,i ( k )) , (7)where Bin ( n, p ) denotes the binomial distribution with parameters n and p . In control-theoreticterms, we model c i ( k ) as the output of a time-varying (due to its dependance on parameters p hx,i ( k )and p hh,i ( k )) stochastic nonlinear observer given by c i ( k ) = g k (∆ s i ( k ) , z i ( k ) , z h,i ( k ) , ξ ) , (8)where ξ captures the randomness in the observation.We now relate the observed number of recoveries ( d i ( k )) with the underlying states in ananalogous manner. Recall that under the SIR epidemic dynamics, the change in the proportion ofrecovered individuals, ∆ r i ( k ) = r i ( k ) − r i ( k −
1) = γ i x i ( k − d i ( k ), which denotes the numberof known removed cases among the known active cases A i ( k − d i ( k ) ∼ Bin ( A i ( k − , γ i ) . (9)In other words, each known active case recovers with probability γ i . When the number of activecases is large, d i ( k ) is approximately equal to γ i A i ( k − The above analysis formally relates the observed quantities with the underlying epidemic states. Ifthe parameters p hx,i ( k ) and p hh,i ( k ) and the number of confirmed cases among the tested high andlow risk populations are known, then maximum likelihood inference of − ∆ s i ( k − τ ) can be computedwithout much difficulty as both p xh,i ( k ) and p xl,i ( k ) are monotonically increasing in − ∆ s i ( k − τ ).Nevertheless, in the rest of the section as well as in our numerical results, we focus on apractically motivated special case where only high risk individuals undergo testing, i.e., z h,i ( k ) = z i ( k ). Note that this has been the practice in many jurisdictions in the world during the ongoingCOVID-19 pandemic [Cohen and Kupferschmidt, 2020]. Following (5), we have c i ( k ) z i ( k ) = (cid:100) p xh,i ( k ) = (cid:34) p hh,i ( k ) p hx,i ( k ) (cid:20) − (cid:99) ∆ s i ( k − τ ) − (cid:21)(cid:35) − = (cid:34) α i ( k ) (cid:20) − (cid:99) ∆ s i ( k − τ ) − (cid:21)(cid:35) − , (10)where α i ( k ) := p hx,i ( k ) p hh,i ( k ) . In other words, we only need to know ratio of the probability of an infectedindividual being high risk (or undergoing testing) and the probability of a healthy individual beinghigh risk (or undergoing testing). The parameter α i ( k ) potentially depends on the number of testscarried out in a given day. Intuition suggests that as z i ( k ) increases, α i ( k ) should converge to 1,while for small values of z i ( k ), it should be very large. The parametric function α i ( k ) = 1+ Fz i ( k ) with F being a large constant is a potential candidate that satisfies the above requirements. Another Modifying the proposed Bayesian approach to incorporate inaccuracy in testing data is beyond the scope of thispaper and will be explored in follow up work. α i ( k ) = N i z i ( k ) where, recall that, N i is the size of sub-population i . Note that if α i ( k ) = 1, i.e., a healthy person is equally likely to undergo testing compared to an infected person,then we have − (cid:99) ∆ s i ( k − τ ) = c i ( k ) z i ( k ) .Given the testing data and assuming that we know α i ( k ), we now discuss how to infer theunderlying states. In particular, for every node i ∈ [ n ], suppose that the daily testing data Ω i ( k )is available to us over the time interval k ∈ [ T + τ, T + τ ]. Accordingly, we define the inferredfraction of new infections at node v i as − (cid:99) ∆ s i ( k ) = (cid:99) ∆ x i ( k ) + (cid:99) ∆ r i ( k ) = (cid:104) − α i ( k + τ ) + α i ( k + τ ) z i ( k + τ ) c i ( k + τ ) (cid:105) − , k ∈ [ T , T ] . (11)However, the above information is not enough to uniquely infer the states. Therefore, we firstassume that for k < T , (cid:98) x i ( k ) = (cid:98) r i ( k ) = 0 and (cid:98) s i ( k ) = 1. In addition, let (cid:99) ∆ r i ( k ) = 0 for k ∈ [ T , T + τ − (cid:99) ∆ x i ( k ) = − (cid:99) ∆ s i ( k ) for k ∈ [ T , T + τ − (cid:98) s i ( k ) = (cid:98) s i ( k −
1) + (cid:99) ∆ s i ( k ) , (cid:98) x i ( k ) = (cid:98) x i ( k −
1) + (cid:99) ∆ x i ( k ) , (cid:98) r i ( k ) = (cid:98) r i ( k −
1) + (cid:99) ∆ r i ( k ) , (12)for k ∈ [ T , T + τ − (cid:99) ∆ r i ( k ) := d i ( k ) A i ( k − (cid:98) x i ( k −
1) = d i ( k ) C i ( k − − D i ( k − (cid:98) x i ( k − , k ∈ [ T + τ, T ] , (13)where (cid:98) x i ( k −
1) is the inferred fraction of infected individuals that is computed recursively following(11) and (12). In the pathological case with A i ( k ) = 0, we assume (cid:99) ∆ r i ( k ) = 0.To summarize, employing (11), (12), and (13) together with the appropriate initialization, wecan infer the underlying epidemic states if testing data is made available over an interval.Note that the inference can be inaccurate when the aforementioned assumptions on the initialconditions are not met. This is indeed the case when testing is carried out after the disease hasalready spread for a while (i.e., when T is large). However, as illustrated in Section 6, simulationsshow that the inferred proportion of infected sub-populations eventually converge to the true values. Remark 2.
The inference problem described above assume knowledge of the parameters α i ( z ) ’s.If these parameters are not known, then one potential approach is to infer the states and solve theleast square problem for different values of α i ( z ) ’s; the one with the smaller cost is likely to be closerto the true underlying α i ( z ) . In the following section, we describe how the inferred data can be used to learn the modelparameters of the SIR epidemic and predict the current state of the epidemic.
Recall from prior discussion that optimal resource allocation to minimize the growth rate of thedisease, or the cost of NPIs, for the SIR epidemic, requires knowledge of the current proportionof susceptible individuals at each node of the network. However, due to the delay τ , testing dataup to time k only suffice to infer the epidemic states up to k − τ . As a result, we need to predict Specifically, we compute (cid:99) ∆ x i ( k ) from our knowledge of (cid:99) ∆ r i ( k ) and − (cid:99) ∆ s i ( k ) and compute (cid:98) x i ( k ) = (cid:98) x i ( k −
1) + (cid:100) ∆ x i ( k ) for k ∈ [ T + τ, T ]. k − τ . Prediction,or forecasting the future trajectory, of the epidemic is also of independent interest (beyond theoptimal resource allocation problem). If the parameters of the epidemic dynamics are known, thenthose can be used to predict the current state of the epidemic from the inferred states up to k − τ .However, at the early stages of the epidemic, there is substantial uncertainty regarding the valuesof β ij and γ i parameters. In this section, we extend the least squares estimation technique proposedin [Chen and Qiu, 2020] for scalar SIR epidemics to learn the unknown β ij and γ i parameters inthe case of networked SIR epidemics.Specifically, suppose at every node i ∈ [ n ], the daily testing data Ω i ( k ) is available to usover the time interval k ∈ [ T + τ, T + τ ]. Following (11)-(13), we can infer the underlyingstates over the interval [ T , T ] from testing data. Subsequently, we can learn the parameters {{ γ i } i ∈ [ n ] , { β ij } ( i,j ) ∈ E } ∈ R n + | E | of the epidemic dynamics by solving the least-square problem:min γ ∈ R n ,β ∈ R | E | n (cid:88) i =1 (cid:104) T (cid:88) k = T (cid:16)(cid:99) ∆ r i ( k ) + (cid:99) ∆ x i ( k ) − (cid:88) j ∈ V i (cid:98) s i ( k ) (cid:98) x j ( k ) β ij (cid:17) + T (cid:88) k = T + τ (cid:16) (cid:99) ∆ r i ( k ) − γ i (cid:98) x i ( k ) (cid:17) (cid:105) . (14)If additional information such as bounds or the sparsity-pattern (specified by the network structure)of β ij parameters are available, those may be incorporated as constraints in (14). We denote theoptimal solutions of (14) as (cid:98) γ i and (cid:98) β ij .The above formulation assumes that the parameters {{ γ i } i ∈ [ n ] , { β ij } ( i,j ) ∈ E } ∈ R n + | E | do notsignificantly change over the interval [ T , T ]. If the parameters change due to NPIs by the author-ities, then the above formulation can be suitably modified to learn the epidemic parameters bothbefore and after the imposition of NPIs by considering suitable sub-intervals during which differentNPIs were in place. Other approaches such as [Chen and Qiu, 2020] define β ij ’s to be parametricfunctions of NPIs. These settings can also be handled by suitably modifying the above formulation. Remark 3.
The formulation in (14) admits a separable structure and as a result, each node v i can learn their respective β ij and γ i values using local information from their own testing data andobtaining the estimates (cid:98) x j ( k ) from their neighbors. The learned parameters of the epidemic dynamics from (14), (cid:98) γ i and (cid:98) β ij , can now be used topredict the current and future state of the epidemic using the inferred states as the initial conditions.In order to avoid introducing additional notation, we use (cid:98) · to also denote the predicted states.Specifically, for k ≥ T , we compute the predicted states as (cid:98) s ( k + 1) = (cid:98) s ( k ) − h diag ( (cid:98) s ( k )) (cid:98) B (cid:98) x ( k ) , (15a) (cid:98) x ( k + 1) = (cid:98) x ( k ) + h diag ( (cid:98) s ( k )) (cid:98) B (cid:98) x ( k ) − h diag ( (cid:98) γ ) (cid:98) x ( k ) , (15b)where (cid:98) B is the matrix with ( i, j )-th entry (cid:98) β ij and diag ( (cid:98) γ ) is a diagonal matrix with diagonal entriesequal to (cid:98) γ i . In the following section, we formulate optimization problems to compute optimal NPIsto control the spread of the epidemic using the predicted state from (15). Recall from Theorem 1 that the growth rate of the infected fraction of the population at time k isgiven by λ A k max (the largest eigenvalue of the matrix I n + h diag ( s ( k )) B − h diag ( γ )). Since A k is afunction of s ( k ), in the absence of perfect knowledge of s ( k ), we use the inferred/predicted value of s ( k ), denoted by (cid:98) s ( k ) (given in (15a)). The goal of the social planner is to control the spread of the11pidemic by choosing the curing rates, i.e., the γ i parameters (which, for instance, correspond todeploying a larger number of healthcare personnel and/or medical equipment) and the contact rates,i.e., the β ij parameters (which correspond to imposing social distancing or lock-down measures).We present two geometric programming formulations to aid the social planner’s decision-makingwith regards to NPIs in a rigorous manner. First, we consider the problem of minimizing theinstantaneous growth rate λ A k max by optimally allocating the NPIs subject to budget constraints.Following analogous arguments in [Han et al., 2015, Preciado et al., 2014], this problem is equivalentto a GP given by: min λ, ¯ γ,w,β λ (16a)s.t. (cid:88) j ∈ V i h (cid:98) s i ( k ) β ij w j /w i + ¯ γ i ≤ λ ∀ i ∈ [ n ] , (16b) (cid:88) ( v i ,v j ) ∈ E f ij ( β ij ) ≤ C , (16c) n (cid:88) i =1 g i (¯ γ i ) ≤ C , (16d) γ l ≤ ¯ γ ≤ γ u , β l ≤ β ≤ β u , (16e) λ ∈ R + , ¯ γ ∈ R n + , w ∈ R n + , β ∈ R | E | + , (16f)where the new variable ¯ γ i := 1 − hγ i is used so that the constraints in (16b) remain posynomials,the constraints in (16e) correspond to bounds on the contact and curing rates, the functions f and g are cost functions of NPIs, the constraints in (16c) and (16d) are budget constraints on NPIswith C and C being the budgets for contact rates and curing rates, respectively.The optimal value λ ∗ k corresponds to the largest eigenvalue of A k with parameters γ ∗ k and β ∗ k ; thelatter denote the optimal curing and contact rates, respectively. The optimal w ∗ k is the eigenvectorcorresponding to λ ∗ k . Furthermore, λ ∗ k is the smallest possible value that can be achieved given thebudget constraints. The above problem is solved repeatedly in an online manner. At time step k + 1, we again obtain (cid:98) s ( k + 1) via feedback, and solve (16) with A k replaced by A k +1 .A problem complementary to (16) is to minimize the cost of NPIs subject to the constraint thatthe growth rate is bounded by λ k . This problem is given by:min ¯ γ,w,β Ψ( β, ¯ γ ) := (cid:88) ( v i ,v j ) ∈ E f ij ( β ij ) + n (cid:88) i =1 g i (¯ γ i ) (17a)s.t. (cid:88) j ∈ V i h (cid:98) s i ( k ) β ij w j /w i + ¯ γ i ≤ λ k ∀ i ∈ [ n ] , (17b) γ l ≤ ¯ γ ≤ γ u , β l ≤ β ≤ β u , (17c)¯ γ ∈ R n + , w ∈ R n + , β ∈ R | E | + . (17d)The above formulation corresponds to imposing NPIs in a cost-optimal manner while ensuring thatthe reproduction number stays below a certain threshold. Remark 4.
Note that in Section 4 we estimated the virus spread parameters while in this sectionwe allow a social planner to control the spread of the virus by setting these parameters. There aretwo ways to interpret this discrepancy: 1) before the social planner is able or decides to exert effortsto mitigate the spread of the virus, they must estimate the state of the system from testing data HITDE ATFR (a) Graph topology analyzed in Section 6. (b) Evolution of infected sub-population
Figure 2: Network topology and evolution SIR epidemic without NPIs.Figure 3: Daily and cumulative synthetic test data for node FR with α = 10. which, given the delay τ , requires estimating the spread parameters, and 2) even after the socialplanner is able and decides to implement preventative measures, the population may not follow therestrictions or they may be less effective (or more extreme) than needed; therefore estimating theactual parameters is necessary. We now evaluate the performance of the proposed framework in the following section via sim-ulations.
In this section, we evaluate the proposed approach on a relatively small network with 5 nodeswith topology given in Fig. 2a. The nodes are roughly modeled after five countries in continentalEurope: France (FR), Germany (DE), Italy (IT), Austria (AT) and Switzerland (CH). Two nodesare neighbors if they share a border. We choose the initially infected proportion to be 0 .
02 at nodeIT and 0 elsewhere. The contact and curing rates are chosen such that Assumption 1 is satisfied.The evolution of the infected sub-population at all the nodes is shown in Figure 2b.
We now generate synthetic test data using the stochastic nonlinear observer proposed in Section 3.We assume that only the high risk (e.g., symptomatic) population undergoes testing, and daily13 igure 4: Daily and cumulative synthetic test data for node FR with α = 1000. number of tests z i ( k ) = 5000 is constant. We generate testing data starting from T = 15 to T = 200 with delay factor τ = 5 by sampling the binomial distribution as shown in (7) and (9).The daily and cumulative numbers of confirmed infections and recoveries as the daily active casesat node FR are shown in Figure 3 and Figure 4 for α FR ( k ) = 10 and α FR ( k ) = 1000, respectively.Recall that c FR ( k ) z FR ( k ) captures the proportion of new infections suitably adjusted by the risk factor α FR ( k ). Recall that α FR ( k ) is the ratio of the probability of an infected individuals being a highrisk and the probability of a healthy individuals being a high risk. As a result, when we only testhigh risk individuals, the proportion of confirmed cases c FR ( k ) z FR ( k ) is much larger when α FR ( k ) is large,and vice versa. In addition, for k ≥ T = 25(i.e., testing is relatively delayed) while all figures on the bottom row have T = 10. The figuresshow that in the absence of early testing, we lack critical information about the initial epidemicstates, and consequently, the inferred values deviate for quite some time from the true values.While the inferred infected proportion is eventually indistinguishable from the true proportion, theinferred recovered proportion remains below the true recovered proportion as the initial informationwas not available, and hence was not incorporated into the inference.The middle and right plots on the top row of Figure 5 differ in the value of α . The figuresshow that as long as α is correctly known and used, the inferred states closely track the truestates. All the figures compare the inferred states and the delayed true states. Next the impactof the parameter estimation and state prediction will be highlighted in the context of optimal NPIcomputation. In order to isolate the performance of the online optimization approaches described in Section 5, wefirst assume that s ( k ) is exactly known, and solve the problems stated in (16) and (17). Following[Preciado et al., 2014], we define the cost functions for NPIs as f ij ( β ij ) = β − ij − β − u,ij β − l,ij − β − u,ij , and g i ( ¯ γ i ) = ¯ γ i − − γ − u,i γ − l,i − γ − u,i ;recall that ¯ γ i = 1 − hγ i . We also set h = 0 .
1. Note that f ij is monotonically decreasing (i.e.,it is costly reduce contact rates) and for β ij ∈ [ β l,ij , β u,ij ], f ij ( β ij ) ∈ [0 , , g i also has similar properties.14 igure 5: Inferred vs. true state trajectories for node IT. Top left and top middle: recovered andinfected sub-population for T = 25, τ = 8 and α = 1000. Top right: infected sub-population for T = 25, τ = 8 and α = 10. Bottom left and middle: T = 10, τ = 8 and α = 1000. Bottom right: T = 10, τ = 1 and α = 10. The upper and lower limits on β ii and β ij are chosen differently to reflect the fact it is easier toreduce contact between individuals from different sub-populations compared to individuals withina sub-population.We first report the results obtained by solving (16) in an online manner for different budgetcombinations. We assume that the initial fraction of infected nodes is 10 − in DE, and 0 for allother nodes. Fig. 6 shows the evolution of the fraction infected and fraction recovered for nodeDE, and the instantaneous growth rate λ ∗ k at the optimal NPIs. The figure shows that higherbudgets lead to “flattening” of the curve of the infected population. The peak of the infectedsub-population is smaller and occurs much later compared to the baseline setting without NPIs.Furthermore, the cumulative number of people who become infected shows significant reductiondepending on the budget. When the budget is sufficiently high, not many people become infectedin the first place. Consequently, the proportion of susceptible individuals remains high, which(perhaps counter-intuitively) results in a larger value of λ ∗ k compared to the case with smallerbudgets for NPIs. Nevertheless, the proportion of infected individuals decreases to 0.The evolution of the infected proportion at node FR and the optimal cost obtained by solving(17) are shown in Fig. 7. We constrain λ A k max to be 1 .
02 for 0 ≤ k ≤ .
95 for k >
100 until the withdrawal of NPIs at k = 200 ,
300 and 500 corresponding to early, intermediateand late withdrawals. When the constraint is 0 .
95, meeting the condition in Theorem 1 for linearconvergence to zero, a significantly larger cost is borne by the social planner.Fig. 7 also shows that early and intermediate withdrawals lead to an exponential increasein the infected proportion, which occurs because, due to the initial interventions, the susceptibleproportion is high and as a result, λ A k max is large, in the absence of NPIs. Therefore, the social plannerincurs a heavy cost for maintaining the growth rate below 1. Unless the NPIs are maintained for a15 igure 6: Evolution of the infected and recovered proportions of node DE and λ ∗ k under optimalbudget-constrained NPIs obtained by solving (16). The baseline scenario refers to the case withoutNPIs.Figure 7: The infected proportion and the optimal cost obtained by solving (17) for differentconstraints on λ A k max . significantly long period to drive down the infected fraction to a negligible amount, there remainsa serious risk of a subsequent flare-up.We are currently investigating the performance of NPIs when the true proportion of the sus-ceptible sub-population is not known, but rather the inferred state information is used to solvethe above optimization problems. Furthermore, inference and parameter estimation using the realtesting data for the countries mentioned above are being carried out. Results pertaining to boththe above factors will be added in a follow up version of this working paper. In this paper, we introduce a closed-loop framework to estimate and control the spread of theSIR epidemic on networks. We first rigorously establish the behavior of the discrete-time SIRepidemic dynamics; specifically that the dynamics is well-defined, the susceptible proportions andthe growth rate decrease monotonically over time, the growth rate eventually falls below 1, and aftersome point the infected proportions of the population converge at least linearly to zero. We alsogeneralize the reproduction number to the network-dependent setting. Furthermore, we incorporateseveral characteristics of real-world testing data in the state estimation task and examine theimpacts of allocating NPIs (such as reducing contact rates and augmenting healthcare equipmentand personnel) by solving suitable GPs in an online manner.16reliminary results reported here provide compelling insights on the behavior of the epidemicdynamics under NPIs. For instance, we show that (i) ignoring the fact that only individuals at ahigher risk of being infected are being tested can have serious impact on inferring the true statesand consequently on the optimal deployment of NPIs, and (ii) premature withdrawal of NPIs couldlead to emergence of a second wave of infections. Nevertheless, this work presents a first attemptat using a closed-loop approach that integrates testing, through forecasting, all the way to control,as opposed to looking at control and forecasting separately. We now summarize several interestingopen problems for the community to explore moving forward. • As mentioned before, we considered the networked SIR model given its strength as a funda-mental and general model. For future work, we would like to expand these ideas to a richerclass of epidemic models with potentially more compartments (such as death, quarantine,exposed, etc) such as the one proposed in [Giordano et al., 2020]. • Similarly, the classical SIR and SEIR models do not fully capture the fact that many individ-uals that become infected with COVID-19 never show any symptoms yet remain infectious.Therefore, a richer class of models with possibly more compartments and/or non-Markoviantransitions should be developed that may better capture the behavior of the COVID-19 pan-demic and other possible outbreaks in the future. • For the estimation part, in Section 4, we have focused on the least-square formulation dueto its simplicity and universality. Other approaches such as Bayesian Markov Chain MonteCarlo [Osthus et al., 2017] may also be used in conjunction with the observation model givenin (7) and (9). • Our numerical evaluations have relied on synthetic data. As was illustrated by the simulations,it is challenging to infer the states when the testing approach is not known and testingdoes not start at the beginning of the outbreak. Ongoing work is focused on inference andlearning from real-data pertaining to COVID-19 pandemic using the proposed framework andproviding insights into the NPIs imposed by various jurisdictions. • Given that the state inference is inaccurate when the deployment of testing is delayed, theimpact of this inaccuracy on the optimal NPI techniques needs to be quantified explicitly inorder to understand the importance of early testing. • Further, if we impose NPIs too soon with the hope of reducing the growth of infections, we maynot have rich enough data to learn the epidemic parameters, which will hinder our forecastingabilities. Without understanding the spread of the disease, it is difficult to know when towithdraw NPIs. Therefore, studying the trade-off between exploration and exploitation isessential. • The success of NPIs depend on the effectiveness with which they are enforced. While the socialplanner computes the optimal contact rates leading to policy guidelines, individuals may ormay not follow them based on their own perceived risk. Therefore, game-theoretic models ofhuman decision-making need to be integrated into the proposed closed-loop framework forimproved prediction and control of the epidemic.Overall, the proposed framework establishes the mathematical foundations for computing optimalNPIs, highlights a number of interesting future research problems and will be a valuable tool forpolicy-makers. This problem was recently investigated for the class of SIS epidemics in [Hota and Sundaram, 2019]. cknowledgement The authors thank Prof. Shreyas Sundaram (Purdue University) and Damir Vrabac (StanfordUniversity) for helpful discussions that contributed to this work. The authors thank Sanket KumarSingh (IIT Kharagpur) for his contributions in data analysis and visualizations.
References
Alberto Abadie, Paolo Bertolotti, Ben Deaner, Arnab Sarker, and Devavrat Shah. Epidemic modeling and estimation,2020. Working Paper, MIT.Alain Barrat, Marc Barthelemy, and Alessandro Vespignani.
Dynamical processes on complex networks . Cambridgeuniversity press, 2008.Stephen Boyd, Seung-Jean Kim, Lieven Vandenberghe, and Arash Hassibi. A tutorial on geometric programming.
Optimization and engineering , 8(1):67, 2007.Francesco Casella. Can the COVID-19 epidemic be managed on the basis of daily data?, 2020. arXiv preprintarXiv:2003.06967.Xiaohui Chen and Ziyi Qiu. Scenario analysis of non-pharmaceutical interventions on global COVID-19 transmissions,2020. arXiv preprint arXiv:2004.04529.Jon Cohen and Kai Kupferschmidt. Countries test tactics in war against COVID-19, 2020.Moez Draief and Laurent Massouli.
Epidemics and rumours in complex networks . Cambridge University Press, 2010.Soheil Eshghi, MHR Khouzani, Saswati Sarkar, and Santosh S Venkatesh. Optimal patching in clustered malwareepidemics.
IEEE/ACM Transactions on Networking , 24(1):283–298, 2014.Giulia Giordano, Franco Blanchini, Raffaele Bruno, Patrizio Colaneri, Alessandro Di Filippo, Angela Di Matteo, andMarta Colaneri. Modelling the COVID-19 epidemic and implementation of population-wide interventions in italy.
Nature Medicine , pages 1–6, 2020.Shuo Han, Victor M Preciado, Cameron Nowzari, and George J Pappas. Data-driven network resource allocation forcontrolling spreading processes.
IEEE Transactions on Network Science and Engineering , 2(4):127–138, 2015.Herbert W Hethcote. The mathematics of infectious diseases.
SIAM Review , 42(4):599–653, 2000.Ashish R Hota and Shreyas Sundaram. Game-theoretic vaccination against networked SIS epidemics and impacts ofhuman decision-making.
IEEE Transactions on Control of Network Systems , 6(4):1461–1472, 2019.Zhiliang Hu, Ci Song, Chuanjun Xu, Guangfu Jin, Yaling Chen, Xin Xu, Hongxia Ma, Wei Chen, Yuan Lin, YishanZheng, et al. Clinical characteristics of 24 asymptomatic infections with COVID-19 screened among close contactsin nanjing, china.
Science China Life Sciences , 63(5):706–711, 2020.Van Sy Mai, Abdella Battou, and Kevin Mills. Distributed algorithm for suppressing epidemic spread in networks.
IEEE Control Systems Letters , 2(3):555–560, 2018.Wenjun Mei, Shadi Mohagheghi, Sandro Zampieri, and Francesco Bullo. On the dynamics of deterministic epidemicpropagation over networks.
Annual Reviews in Control , 44:116–128, 2017.Joel C Miller. A note on the derivation of epidemic final sizes.
Bulletin of mathematical biology , 74(9):2125–2141,2012.Cameron Nowzari, Victor M Preciado, and George J Pappas. Analysis and control of epidemics: A survey of spreadingprocesses on complex networks.
IEEE Control Systems , 36(1):26–46, 2016.Masaki Ogura and Victor M Preciado. Efficient containment of exact SIR Markovian processes on networks. In , pages 967–972. IEEE, 2016. asaki Ogura and Victor M Preciado. Optimal containment of epidemics in temporal and adaptive networks. In Temporal Network Epidemiology , pages 241–266. Springer, 2017.Masaki Ogura, Victor M Preciado, and Naoki Masuda. Optimal containment of epidemics over temporal activity-driven networks.
SIAM Journal on Applied Mathematics , 79(3):986–1006, 2019.Dave Osthus, Kyle S Hickmann, Petrut¸a C Caragea, Dave Higdon, and Sara Y Del Valle. Forecasting seasonalinfluenza with a state-space SIR model.
The Annals of Applied Statistics , 11(1):202, 2017.Miyo Ota. Will we see protection or reinfection in COVID-19?, 2020.Philip E Par´e, Ji Liu, Carolyn L Beck, Barrett E Kirwan, and Tamer Ba¸sar. Analysis, estimation, and validation ofdiscrete-time epidemic processes.
IEEE Transactions on Control Systems Technology , 28(1):79–93, 2020.Romualdo Pastor-Satorras, Claudio Castellano, Piet Van Mieghem, and Alessandro Vespignani. Epidemic processesin complex networks.
Reviews of Modern Physics , 87(3):925, 2015.Victor M Preciado, Michael Zargham, Chinwendu Enyioha, Ali Jadbabaie, and George J Pappas. Optimal resourceallocation for network protection against spreading processes.
IEEE Transactions on Control of Network Systems ,1(1):99–108, 2014.Eduardo Ram´ırez-Llanos and Sonia Mart´ınez. Distributed discrete-time optimization algorithms with applicationsto resource allocation in epidemics control.
Optimal Control Applications and Methods , 39(1):160–180, 2018.Peter X Song, Lili Wang, Yiwang Zhou, Jie He, Bin Zhu, Fei Wang, Lu Tang, and Marisa Eisenberg. An epidemio-logical forecast model and software assessing interventions on COVID-19 epidemic in china, 2020. medRxiv.Richard S. Varga.
Matrix Iterative Analysis . Springer-Verlag, 2000.Yan Wan, Sandip Roy, and Ali Saberi. Network design problems for controlling virus spread. In , pages 3925–3932. IEEE, 2007.Yan Wan, Sandip Roy, and Ali Saberi. Designing spatially heterogeneous strategies for control of virus spread.
IETSystems Biology , 2(4):184–201, 2008.Nicholas J Watkins, Cameron Nowzari, and George J Pappas. Robust economic model predictive control of continuous-time epidemic processes.
IEEE Transactions on Automatic Control , 65(3):1116–1131, 2019.Gul Zaman, Yong Han Kang, and Il Hyo Jung. Stability analysis and optimal vaccination of an SIR epidemic model.
BioSystems , 93(3):240–249, 2008.
A Appendix: Proof of Theorem 1.
Proof.
We present the proof for each part of the theorem, starting with 1).1) We prove this result by induction. By assumption s i (0) , x i (0) , r i (0) ∈ [0 ,
1] and s i (0)+ x i (0)+ r i (0) = 1 for all i ∈ [ n ]. Therefore, by Assumption 1 and (3a), s i (1) ≥ (1 − h (cid:80) nj =1 β ij ) s i (0) ≥
0. Since h (cid:2) − s i (0) (cid:80) nj =1 β ij x j (0) (cid:3) ≤ s (1) ≤ s (0) ≤
1. By Assumption 1 and (3b), x i (1) ≥ (1 − hγ i ) x i (0) ≥
0. Since x j (0) ≤ j ∈ [ n ], by (3b) and Assumption 1, x i (1) ≤ x i (0) + s i (0) h (cid:80) nj =1 β ij ≤ x i (0)+ s i (0) ≤
1. By (3c) and the non-negativity of h and γ i , and since x i (0) ≥ r i (1) ≥ r i (0) ≥
0. By (3b) and Assumption 1, we have r i (1) ≤ r i (0) + x i (0) ≤
1. Addingup (3a)-(3c), gives that s i (1) + x i (1) + r i (1) = s i (0) + x i (0) + r i (0), which by assumption equals 1.Now assume for an arbitrary k , s i ( k ) , x i ( k ) , r i ( k ) ∈ [0 ,
1] and s i ( k )+ x i ( k )+ r i ( k ) = 1. Followingthe exact same arguments as for the base case except replacing 0 with k and 1 with k + 1, it can beshown that s i ( k + 1) , x i ( k + 1) , r i ( k + 1) ∈ [0 ,
1] and s i ( k + 1) + x i ( k + 1) + r i ( k + 1) = 1. Therefore,by induction, s i ( k ) , x i ( k ) , r i ( k ) ∈ [0 ,
1] and s i ( k ) + x i ( k ) + r i ( k ) = 1 for all k ≥ i ∈ [ n ].2) By 1) and the non-negativity of the β ij ’s, it follows that h (cid:2) − s i ( k ) (cid:80) nj =1 β ij x j ( k ) (cid:3) ≤ k ≥
0. Therefore, from (3a), we have s i ( k + 1) ≤ s i ( k ).19) Since the rate of change of s ( k ), − h diag ( s ( k )) Bx ( k ), is non-positive for all k ≥ s ( k )is lower bounded by zero, by 1), we conclude that lim k →∞ s ( k ) exists. Thereforelim k →∞ − h diag ( s ( k )) Bx ( k ) = n , (18)where n is the vector of dimension n with all entries equal to 0. Therefore, lim k →∞ x ( k +1) − x ( k ) =lim k →∞ − hγx ( k ). Thus, by the assumption that hγ i > i ∈ [ n ], lim k →∞ x i ( k ) = 0 for i ∈ [ n ].4) Recall that by Assumption 1, A k is an irreducible non-negative matrix and thus by thePerron-Frobenius Theorem for irreducible non-negative matrices [Varga, 2000, Theorem 2.7 andLemma 2.4], λ A k max = ρ ( A k ), where ρ ( · ) indicates the spectral radius. By [Varga, 2000, Theorem 2.7], ρ ( A k ) increases when any entry of A k increases. Therefore, by 2) and since A k is defined as I n + h diag ( s k ) B − h diag ( γ ), we have that ρ ( A k ) ≥ ρ ( A k +1 ) , which implies λ A k max ≥ λ A k +1 max .
5) There are two possible types of equilibria for the SIR model: i) lim k →∞ s ( k ) = n , orii) lim k →∞ s ( k ) = s ∗ (cid:54) = n . We explore the two cases separately.i) If lim k →∞ s ( k ) = n , then the rate of change of x ( k ) converges to − hγx ( k ). Therefore, bythe definition of λ A k max , there exists a ¯ k such that λ A k max < k ≥ ¯ k .ii) If lim k →∞ s ( k ) = s ∗ (cid:54) = n , then, by 3), for any ( s (0) , x (0) , r (0)) the system dynamics convergeto some equilibria of the form ( s ∗ , n , n − s ∗ ), for some nonzero s ∗ . Define (cid:15) s ( k ) := s ( k ) − s ∗ and (cid:15) x ( k ) := x ( k ) − n . (19)By 2) and 1), respectively, we know that (cid:15) s ( k ) ≥ (cid:15) x ( k ) ≥ k ≥
0. Furthermore, weknow that (cid:15) s ( k + 1) ≤ (cid:15) s ( k ) for all k ≥
0, lim k →∞ (cid:15) s ( k ) = n , and lim k →∞ (cid:15) x ( k ) = n .Linearizing the dynamics of (cid:15) s ( k ) and (cid:15) x ( k ) around ( s ∗ , n ) gives (cid:15) s ( k + 1) = (cid:15) s ( k ) − h diag ( s ∗ ) B(cid:15) x ( k ) , (20a) (cid:15) x ( k + 1) = (cid:15) x ( k ) + h diag ( s ∗ ) B(cid:15) x ( k ) − h diag ( γ ) (cid:15) x ( k ) . (20b)Let λ A ∗ max be the maximum eigenvalue of I n + h diag ( s ∗ ) B − h diag ( γ ) with corresponding normalizedleft eigenvector w ∗ , that is, w ∗(cid:62) ( I n + h diag ( s ∗ ) B − h diag ( γ )) = λ A ∗ max w ∗(cid:62) . (21)If λ A ∗ max >
1, then the system in (20) is unstable. Therefore, by Lyapunov’s Indirect Method,lim k →∞ ( (cid:15) s ( k ) , (cid:15) x ( k )) (cid:54) = ( s ∗ , n ), which is a contradiction.Now consider the case where λ A ∗ max = 1. Left multiplying the equation for x ( k + 1) in (4b) by w ∗(cid:62) and using (19) and (21) gives w ∗(cid:62) (cid:15) x ( k + 1) = w ∗(cid:62) [ I n + h diag ( (cid:15) s ( k ) + s ∗ ) B − h diag ( γ )] (cid:15) x ( k )= λ A ∗ max w ∗(cid:62) (cid:15) x ( k ) + w ∗(cid:62) h diag ( (cid:15) s ( k )) B(cid:15) x ( k )= w ∗(cid:62) (cid:15) x ( k ) + w ∗(cid:62) h diag ( (cid:15) s ( k )) B(cid:15) x ( k ) , where the last equality holds since λ A ∗ max = 1. Thus, w ∗(cid:62) ( (cid:15) x ( k + 1) − (cid:15) x ( k )) = w ∗(cid:62) h diag ( (cid:15) s ( k )) B(cid:15) x ( k ) ≥ , k →∞ x ( k ) = n . Therefore, there exists a ¯ k such that λ A k max < k ≥ ¯ k .6) Since, by 5), there exists a ¯ k such that λ A k max < k ≥ ¯ k , and we know that λ A k max = ρ ( A k ) ≥ k →∞ (cid:107) x ( k + 1) (cid:107)(cid:107) x ( k ) (cid:107) = (cid:107) A k x ( k ) (cid:107)(cid:107) x ( k ) (cid:107) = λ A k max < . (22)Therefore, for k ≥ ¯ k , x ( k ) converges linearly to n .Further, if lim k →∞ s i ( k ) = 0 and hγ i = 1 for all i ∈ [ n ], there exists ¯¯ k , such that λ A k max = 0.Therefore, by (22), x i ( k ) converges superlinearly to 0 for all k ≥ ¯¯ kk