Estimation and Distributed Eradication of SIR Epidemics on Networks
Ciyuan Zhang, Humphrey Leung, Brooks Butler, Philip. E. Paré
11 Estimation and Distributed Eradication of SIR Epidemics on Networks
Ciyuan Zhang, Humphrey Leung, Brooks Butler, and Philip. E. Paré*
Abstract — This work examines the discrete-time networkedSIR (susceptible-infected-recovered) epidemic model, where theinfection and recovery parameters may be time-varying. Weprovide a sufficient condition for the SIR model to converge tothe set of healthy states exponentially. We propose a stochasticframework to estimate the system states from observed testingdata and provide an analytic expression for the error of theestimation algorithm. Employing the estimated and the truesystem states, we provide two novel eradication strategies thatguarantee at least exponential convergence to the set of healthystates. We illustrate the results via simulations over northernIndiana, USA.
I. I
NTRODUCTION
As of February 2021, the COVID-19 virus has claimed2.4 million lives and infected 110 million individuals world-wide [1]. Lack of effective treatments, high contagionrates [2], long incubation periods [3]–[6], and asymptomaticcases [7]–[10] pose significant challenges in containing anderadicating pandemics. Recent pandemics, including gonor-rhea [11], Ebola [12], and COVID-19 [13], have acceleratedthe development of infection models. The main goal ofepidemic model development is to identify conditions toeradicate the pathogen, and leverage the knowledge of theseconditions to design mitigation strategies. Various infec-tion models have been proposed, based on characteristicsof individual pathogens, and studied in the literature, in-cluding susceptible-infected-susceptible (SIS), susceptible-infected-removed (SIR), and susceptible-infected-removed-susceptible (SIRS) [14], [15]. In this paper, we focus on theSIR epidemic model. We aim to expand on the SIR model,by exploring mutating viruses over networks, estimation ofthe underlying states, and distributed eradication strategies.The patchwork response to COVID-19 [16] gives rise tosusceptible community subpopulations, with heterogeneoustime-varying factors not previously explored by the SIRmodel. Extensions on the SIS model, studied in [17]–[19], augment the compartmental epidemic models origi-nated in [20] to include interactions between subpopulationsof susceptible communities. Additionally, various advancedepidemic models consider time-varying factors [21]–[28]. Inthis paper, we establish sufficient conditions for the set ofhealthy states of a networked time-varying SIR model to beglobally exponentially stable (GES). These equilibrium statesare not unique, as the final susceptible and removed states *Ciyuan Zhang, Humphrey Leung, Brooks Butler, and Philip. E. Paré arewith the School of Electrical and Computer Engineering at Purdue Univer-sity. Emails: {zhan3375, leung61, brooksbutler, philpare}@purdue.edu.*This work was funded in part by the C3.ai Digital TransformationInstitute sponsored by C3.ai Inc. and the Microsoft Corporation and in partby the National Science Foundation, grants NSF-CNS are dependent on the initial conditions and the time-varyinginfectious and healing parameters.The delay in onset of COVID-19 symptoms [3]–[6], largeasymptomatic populations estimated between − [7]–[10], and delay in test results [29] compromise the ability foraccurate estimation of current infection states. An estimationalgorithm that incorporated a constant delay between thechange in infection proportion and testing data was intro-duced in [30]. They studied the inference problem by usinga Bayesian approach. Inspired by the delay characterizationsuggested in [30], we propose a stochastic delay to modelthe unpredictability of the COVID-19 virus and testingstrategies. We have developed methods for estimating theunderlying epidemic states from testing data with a delaysampled from a geometric distribution, which cannot becompletely filtered by the method suggested in [30]. Thegeometric delay model accounts for the stochastic effect ofindividuals failing to get tested immediately after exposure.We study the aggregated effect of each individual delay onthe trajectory of confirmed cases and devise a method forestimating the underlying epidemic states of an SIR modelfrom these delayed measurements. We also investigate theproposed method’s estimation error, which provides insightsfor achieving an accurate estimation of the system states. Wethen employ this more realistic estimation to strategicallyeradicate a disease.As proven in [30], the SIR epidemic model convergesto a healthy state, however, an exponential convergence isnot shown. Combining the modeling and inference approachallows us to develop two distributed control strategies ca-pable of eradicating epidemic spread exponentially, at anequilibrium with a higher proportion of susceptible popu-lation. Decreasing the removed (recovered) and increasingthe susceptibility proportion is of exceptional importance forthe COVID-19 pandemic, as long term and severe healthcomplications have been documented in the recovered pop-ulations, including impaired cognition [31] and damage tocardiac tissue [32]. Our main result shows that by applyingthe estimated and the true susceptible states of each node,our proposed eradication strategies will guarantee globalexponential stability of a healthy state of the overall network. A. Paper Contributions
We summarize the main contributions of this paper asfollows: • We establish sufficient conditions for global exponentialstability of the set of healthy states; see Theorem 1. • We propose a stochastic framework which estimatesthe trajectories of the system states of the networked a r X i v : . [ ee ss . S Y ] F e b SIR model from testing data. Furthermore, we provideanalytical expressions for the error of the estimationalgorithm we propose; see Prop. 1. • We propose two distributed eradication strategies foradjusting healing rates, one that is based on the truesystem states and the other one is based on the inferredsystem states. Both methods guarantee that the virus iseradicated within exponential time; see Theorem 2 andCorollary 3.
B. Paper Outline
We organize this paper as follows: Section II lays downsome basic assumptions and restates the well-known SIRmodel in the networked fashion, and it presents the mainproblems studied in this paper. Section III first recalls pre-liminary results that are essential for stability analysis andthen discusses the sufficient conditions for global exponentialstability of a healthy state of the networked time-varyingSIR models. Section IV covers the proposed techniques ofestimating hidden epidemics states with the stochastic delayof tested individuals and testing data. Section V coversthe two distributed control strategies which ensure that thesystem converges to a healthy state in at least exponentialtime. Section VI illustrates the results from Section IV andV with numerical simulations. Finally, in Section VII, wesummarize the main conclusions of this paper and discusspotential future directions.
C. Notation
We denote the set of real numbers, the non-negative inte-gers, and the positive integers as R , Z ≥ , and Z ≥ , respec-tively. For any positive integer n , we have [ n ] = { , , ..., n } .The spectral radius of a matrix A ∈ R n × n is ρ ( A ) . Adiagonal matrix is denoted as diag ( · ) . The transpose of avector x ∈ R n is x (cid:62) . The Euclidean norm is denoted by (cid:107)·(cid:107) .We use I to denote identity matrix. We use and to denotethe vectors whose entries all equal 0 and 1, respectively. Thedimensions of the vectors are determined by context. Givena matrix A , A (cid:31) (resp. A (cid:23) ) indicates that A is positivedefinite (resp. positive semidefinite), whereas A ≺ (resp. A (cid:22) ) indicates that A is negative definite (resp. negativesemidefinite). Let G = ( V , E ) denote a graph or networkwhere V = { v , v , ..., v n } is the set of subpopulations, and E ⊆ V × V is the set of edges. We denote the expectationof a random variable as IE[ · ] .II. M ODEL AND P ROBLEM F ORMULATION
Consider a time-varying epidemic network of n subpopu-lations, where the size of subpopulation v i is N i ∈ Z > ,and the infection rates and healing rates could be time-varying. We denote β ij ( t ) as the infection rate from node v j to node v i at time t , we denote γ i ( t ) as the healing rateof node v i at time t . The proportions of the subpopulationat node v i which are susceptible, infected, and recovered attime t are denoted by s i ( t ) , x i ( t ) and r i ( t ) , respectively. The deterministic continuous-time evolution of the SIR epidemicis given by ˙ s i ( t ) = − s i ( t )[ n (cid:88) j =1 β ij ( t ) x j ( t )] , (1a) ˙ x i ( t ) = s i ( t )[ n (cid:88) j =1 β ij ( t ) x j ( t )] − γ i ( t ) x i ( t ) , (1b) ˙ r i ( t ) = γ i ( t ) x i ( t ) , ∀ i ∈ [ n ] . (1c)We now state the discrete-time SIR epidemic dynamicsobtained through Euler discretization of (1). For a smallsampling time h > , the discrete-time evolution of the SIRepidemic is given by s i ( k + 1) = s i ( k ) + h [ − s i ( k ) n (cid:88) j =1 β ij ( k ) x j ( k )] , (2a) x i ( k + 1) = x i ( k ) + h [ s i ( k ) n (cid:88) j =1 β ij ( k ) x j ( k ) − γ i ( k ) x i ( k )] , (2b) r i ( k + 1) = r i ( k ) + hγ i ( k ) x i ( k ) . (2c)Equation (2b) can be rewritten as x ( k + 1) = x ( k ) + h [ S ( k ) B ( k ) − Γ( k )] x ( k ) , (3)where S ( k ) = diag ( s ( k )) , B ( k ) is the matrix with ( i, j ) thentry β ij ( k ) , and Γ( k ) = diag ( γ i ( k )) . The spread of a virusover a network can be captured using a graph G = ( V , E ) ,where E = { ( v i , v j ) | β ij ( k ) (cid:54) = 0 } is the set of directed edges.We make the following assumptions in order for thesystem in (2) to be well defined. Assumption 1.
For every i ∈ [ n ] , hγ i ( k ) > and ∀ j ∈ [ n ] , β ij ( k ) ≥ , for every k ∈ Z ≥ . Assumption 2.
For every i ∈ [ n ] , hγ i ( k ) ≤ and h (cid:80) j β ij ( k ) ≤ , for every k ∈ Z ≥ . We have the following result which shares the same ideaas the time-invariant model, proven in [30].
Lemma 1.
Suppose s i (0) , x i (0) , r i (0) ∈ [0 , , s i (0) + x i (0) + r i (0) = 1 , and Assumptions 1 and 2 hold. Then,for all k ∈ Z ≥ ,1) s i ( k ) , x i ( k ) , r i ( k ) ∈ [0 , ,2) s i ( k ) + x i ( k ) + r i ( k ) = 1 , and3) s i ( k + 1) ≤ s i ( k ) . Definition 1.
We define the set of healthy states of (2) as { s ∗ i ( k ) , x ∗ i ( k ) , r ∗ i ( k ) : i ∈ [ n ] , k ∈ Z ≥ } , where x ∗ i ( k ) = 0 , s ∗ i ( k ) ∈ [0 , , and r ∗ i ( k ) ∈ [0 , for all i ∈ [ n ] . Given a network that is infected by a virus, our goal is toguarantee that each subpopulation v i converges to the set ofhealthy states in exponential time regardless of the initialconditions of the each subpopulations. We now officiallystate the questions being studied in this paper:(i) For the system with dynamics given in (3), under whatcondition is the set of healthy states, i.e., x ( k ) = ,global exponentially stable (GES)? (ii) Given the testing data, how can the stochastic frame-work be constructed to estimate the susceptible, in-fected, and recovered proportions, denoted by (cid:98) s i ( k ) , (cid:98) x i ( k ) , and (cid:98) r i ( k ) , respectively, for each subpopulation v i in the network?(iii) What is the estimation error of the stochastic frameworkthat we proposed?(iv) Given the knowledge of the conditions that ensure GESof a healthy state, i.e., x ( k ) = and (cid:98) s i ( k ) inferredfrom testing data, how can we devise dynamic controlalgorithms which apply new healing rates (cid:98) γ i ( k ) to eachagent in (3) so that the epidemic is eradicated with afaster rate of convergence than the rate of exponential?III. S TABILITY A NALYSIS
This section presents conditions that ensure global ex-ponential stability of the set of healthy states. First, weintroduce some preliminaries and then we present our mainanalysis results.
A. Preliminaries
In this subsection, we recall results that are crucial forunderstanding the rest of the paper.
Lemma 2. [33] Suppose that M is a nonnegative matrixwhich satisfies ρ ( M ) < . Then there exists a diagonalmatrix P (cid:31) such that M (cid:62) P M − P ≺ . Consider a system described as follows: x ( k + 1) = f ( k, x ( k )) . (4) Definition 2.
An equilibrium point of (4) is GES is thereexist positive constants α and ω , with ≤ ω < , such that (cid:107) x ( k ) (cid:107) ≤ α (cid:107) x ( k ) (cid:107) ω ( k − k ) , ∀ k, k ≥ , ∀ x ( k ) ∈ R n . (5)We recall a sufficient condition for GES of an equilibriumof (4) from [34]. Lemma 3. [34, Theorem 28] Suppose there exists a function V : Z + × R n → R , and constants a, b, c > and p > suchthat a (cid:107) x (cid:107) p ≤ V ( k, x ) ≤ b (cid:107) x (cid:107) p , ∆ V ( k, x ) := V ( x ( k +1)) − V ( x ( k )) ≤ − c (cid:107) x (cid:107) p , ∀ k ∈ Z ≥ , and ∀ x ( k ) ∈ R n ,then x ( k ) = is a globally exponential stable equilibriumof (4) . Lemma 4. [35, Theorem 23.3] Under the assumption ofLemma 3, the rate of convergence to the origin is upperbounded by an exponential rate of (cid:112) − ( c/b ) ∈ [0 , ,where b and c are defined in Lemma 3. Note that a healthy state of the system in (2) is GES ifAssumptions 1 and 2 hold and the condition in Definition 2and Lemma 3 are satisfied for all x ( k ) ∈ [0 , n , since thisis the domain where the model is well defined. B. Global Exponential Stability of the Healthy States
In this subsection, we present sufficient conditions for theglobal exponential stability of the set of healthy states ofthe system. We find the conditions by analyzing the spectralradius of the state transition matrix of (2b). We define M ( k ) = I − h Γ( k ) + hB ( k ) , (6) ˆ M ( k ) = I + h [ S ( k ) B ( k ) − Γ( k )] . (7)Notice that ˆ M ( k ) is the state transition matrix of (2b) andit can be written that ˆ M ( k ) = M ( k ) − h ( I − S ( k )) B ( k ) . (8)We use M ( k ) and (8) to illustrate the sufficient conditionsfor the GES of the set of healthy states in the subsequenttheorem. Theorem 1.
Given Assumptions 1 and 2, suppose for all k ∈ Z ≥ , B ( k ) is symmetric. If sup k ∈ Z ≥ ρ ( M ( k )) < ,then the set of healthy states of (2) is GES.Proof. See Appendix.Recall from the previous result that M ( k ) is a nonnegativematrix which satisfies sup k ∈ Z ≥ ρ ( M ( k )) < , such that M (cid:62) ( k ) Q ( k +1) M ( k ) − Q ( k ) ≺ , where Q ( k ) is a diagonalmatrix defined in the Lyapunov function: V ( k, x ) = x (cid:62) Q ( k ) x. (9) Corollary 1.
Under the assumptions of Theorem 1,the rate of convergence to a healthy state is upperbounded by an exponential rate of (cid:113) − σ σ , where σ = max k ∈ Z ≥ λ max ( Q ( k )) , σ = max k ∈ Z ≥ λ min [ Q ( k ) − M ( k ) (cid:62) Q ( k + 1) M ( k )] .Proof. See Appendix.
Remark 1.
Notice that in Theorem 1, sup k ∈ Z ≥ ρ ( M ( k )) < is the key condition which ensures the set of healthystates of (2) is GES. We can interprete ρ ( M ) in the contextof epidemiology as the basic reproduction number of thevirus over the network. Particularly, Theorem 1 affirms thatgiven the time-varying parameters of the network satisfythe condition provided, the mutating virus will exponentiallyconverge to the set of healthy states. In this section, we found the conditions that ensure ex-ponential convergence to the set of healthy states of (2b)in Theorem 1. This answers question (i) in Section II. Thestability condition can help the policymakers reallocate themedical resources, staff etc. which leads to modifying theparameters in (2) so that the spreading of the virus stopscompletely. One of the other factors that will assist indecision making is the COVID-19 observed testing data.IV. S
TATE E STIMATION FROM T ESTING D ATA
In this section, we study how to estimate the epidemicstates ( s ( k ) , x ( k ) , r ( k )) from testing data in order to designa feedback controller in the following section. One of thechallenges of estimating the underlying system states is that Ω i ( k ) = (cid:0) C i ( k ) , D i ( k ) (cid:1) → (cid:0) c i ( k ) , d i ( k ) (cid:1)(cid:98) Θ i ( k − τ i ) = (cid:0) (cid:98) s i ( k − τ i ) , (cid:98) x i ( k − τ i ) , (cid:98) r i ( k − τ i ) (cid:1)(cid:98) Θ i ( k ) = (cid:0) (cid:98) s i ( k ) , (cid:98) x i ( k ) , (cid:98) r i ( k ) (cid:1) FIGURE
1: Estimation of System States from Testing Datathe testing data on a given day does not capture the newinfections on the same day. Instead, the testing data is adelayed representation of the change in the system. Char-acterizing the delay of each individual is difficult, becausethe delay is determined by numerous factors such as theincubation period of COVID-19, the duration of obtainingtest results, the willingness of each individual to get tested,etc. Therefore, we propose a stochastic framework in thissection to capture the factors which cause the testing delay.
Definition 3.
The testing delay τ i is the length of timebetween when an individual from subpopulation v i is infectedand when their positive test result is reported. In our discrete-time model, we assume that τ i ∈ Z ≥ .We model the testing delay of each infected individual τ i astwo aggregate components to represent the uncertainty in thetesting process: τ i = η i + Y i , (10)where η i ∈ Z ≥ is a constant and Y i is sampled froma discrete time random variable whose measurable spaceis Z ≥ . Remark 2. In (10) , the constant component η i can beinterpreted as the length of time needed to acquire testingresults. The random variable can be interpreted as theincubation period and/or the amount of time that it takesan individual to get tested after becoming infected. First, we denote the set of estimated system states forsubpopulation v i at time k as (cid:98) Θ i ( k ) = (cid:0) (cid:98) s i ( k ) , (cid:98) x i ( k ) , (cid:98) r i ( k ) (cid:1) ,we denote the set of testing data recorded at time k tobe Ω i ( k ) = (cid:0) C i ( k ) , D i ( k ) (cid:1) , where C i ( k ) is the numberof confirmed cases at time k , and D i ( k ) represents thenumber of removed (recovered) cases at time k . In addition,the cumulative number of confirmed and removed cases atnode v i are written as C i ( k ) = (cid:80) kj =0 C i ( j ) and D i ( k ) = (cid:80) kj =0 D i ( j ) , respectively. Therefore, the number of activecases is calculated by A i ( k ) = C i ( k ) − D i ( k ) . Recall thatthe size of each subpopulation is N i ; we define c i ( k ) = C i ( k ) N i and d i ( k ) = D i ( k ) N i as the proportion of daily confirmed casesand removal, respectively. Note that c i ( k ) , d i ( k ) ∈ [0 , , for all i ∈ [ n ] , k ∈ Z ≥ . The estimation procedure is illustratedin Fig. 1.We then study how to relate c i ( k ) to the underlyingstates. We define a vector space Π T as the space of all theproportions of daily number of confirmed cases from timestep k = T to time step k = T + 1 . We define Ξ T as thevector of all the decreases in the proportion of susceptibleindividuals, − ∆ s i ( k ) , from time step k = T to time step k = T + 1 . We denote Φ( T , T ) as the transfer matrixwhich results in Π T = Φ( T , T )Ξ T , (11)where Φ( T , T ) is a ( T − T + 2) × ( T − T + 2) matrix,which depends on the SIR dynamics, the testing strategies,and the delay.When the testing delay is a constant, i.e., τ i = η i > ,the only non-zero entries in Φ( T , T ) are: Φ l + η i ,l = 1 , l ∈ [ T − T + 2 − η i ] . Since for all k ∈ [ T + η i , T + 1] , wecan write c i ( k ) as c i ( k ) = − ∆ s i ( k − η i ) . (12)When the delay η i = 0 , the transfer matrix Φ( T , T ) = I .Because for all k ∈ [ T , T + 1] , we can write that c i ( k ) = − ∆ s i ( k ) . (13)We now propose a stochastic testing framework to capturethe delay between when an individual is infected and whenthey receive a positive test result. We first let η i = 0 in (10),without the loss of generality. Furthermore, we assume thateach infected individual at node v i has an equal probability p xi ∈ (0 , of receiving a diagnostic test each day startingfrom the day after they are infected. Therefore, we model Y i in (10) as a random variable following the geometricdistribution, with the probability of an infected individualacquiring a positive test δ days after infection being: P ( Y i = δ ) = p xi (1 − p xi ) δ − (14)for δ ∈ Z ≥ . The geometric distribution of the testing delaymodels the number of days before an infected individual ob-tains a diagnostic test which represents the incubation periodof COVID-19 and/or the unwillingness of each individualgetting a test. We assume that the delay of each infectedindividual’s distribution is i.i.d. (independent and identicallydistributed) from others. Furthermore, we assume that aninfected individual can be tested only once. Based on [36]and [37], we assume that even if an individual recovers fromCOVID-19, their antibody tests will still give positive results.We also assume that all the tests generate accurate results.We now relate the proportion of confirmed cases c i ( k ) with the underlying states of the system. We define a binaryrandom variable X i ( ν ) with X i ( ν ) = 1 (resp. X i ( ν ) = 0 ) ifa randomly chosen individual from subpopulation v i becameinfected at time ν . It can be written that X i ( ν ) = (cid:40) w.p. − ∆ s i ( ν )0 w.p. s i ( ν ) , (15)where, from (2), − ∆ s i ( ν ) = s i ( ν − − s i ( ν ) = hs i ( ν − (cid:80) j β ij x j ( ν − ≥ for all ν ≥ . We define the binary random variable T i ( µ, δ ) , with T i ( µ, δ ) = 1 if a randomly chosen individual acquired apositive test at time µ and was infected δ days before µ . Nowwe rewrite ν , in (15), as µ − δ . From (14), the conditionalprobability P ( T i ( µ, δ ) = 1 |X i ( µ − δ ) = 1) is given by thegeometric probability mass function (pmf) p xi (1 − p xi ) δ − andrepresents the probability of an infected individual acquiringa positive test specifically δ days after being infected. Hence,the joint pmf of the two random variables X i ( µ − δ ) , T i ( µ, δ ) is written as: P X i , T i ( µ − δ, µ ) = P ( X i ( µ − δ ) ∩ T i ( µ, δ )) , (16)which is interpreted as the probability that a randomly chosenindividual became infected at time µ − δ and acquired apositive test at time µ , where µ − δ, µ ∈ [ T , T ] .Therefore, the joint pmf P X i , T i ( µ − δ, µ ) is calculated as: P X i , T i ( X i ( µ − δ ) = 1 ∩ T i ( µ, δ ) = 1)= P ( T i ( µ, δ ) = 1 |X i ( µ − δ ) = 1) P ( X i ( µ − δ ) = 1)= p xi (1 − p xi ) δ − [ − ∆ s i ( µ − δ )] , (17) P X i , T i ( X i ( µ − δ ) = 0 ∩ T i ( µ, δ ) = 1)= P ( T i ( µ, δ ) = 1 |X i ( µ − δ ) = 0) P ( X i ( µ − δ ) = 0)= 0[1 + ∆ s i ( µ − δ )] = 0 , (18)since we assume that the test results are accurate. Similarly, P X i , T i ( X i ( µ − δ ) = 1 ∩ T i ( µ, δ ) = 0)= P ( T i ( µ, δ ) = 0 |X i ( µ − δ ) = 1) P ( X i ( µ − δ ) = 1)= [1 − p xi (1 − p xi ) δ − ][ − ∆ s i ( µ − δ )] ,P X i , T i ( X i ( µ − δ ) = 0 ∩ T i ( µ, δ ) = 0)= P ( T i ( µ, δ ) = 0 |X i ( µ − δ ) = 0) P ( X i ( µ − δ ) = 0)= 1 + ∆ s i ( µ − δ ) . Let W i ( µ ) be the marginal distribution of T i ( µ, δ ) over theset of feasible delays, δ , with its pmf being the probabilityof a random individual acquiring a positive test at time µ : P W i ( W i ( µ ) = 1)= µ − T (cid:88) δ =1 P X i , T i ( X i ( µ − δ ) ∩ T i ( µ, δ ) = 1)= µ − T (cid:88) δ =1 p xi (1 − p xi ) δ − [ − ∆ s i ( µ − δ )] , by combining (17) and (18). Therefore, the number ofconfirmed cases at time k is calculated as C i ( k ) = IE[ N i (cid:88) l =1 W i ( k )]= N i (cid:88) l =1 IE[ W i ( k )] (19) = N i k − T (cid:88) δ =1 p xi (1 − p xi ) δ − [ − ∆ s i ( k − δ )] , (20) where (19) holds because of the linearity of expectationand since the testing delays are i.i.d. Hence, by combining c i ( k ) = C i ( k ) N i and (20) for all k ∈ [ T + 1 , T + 1] , thetransfer matrix Φ( T , T ) in (11) is written as Φ( T , T )= . . . p xi . . . p xi (1 − p xi ) p xi . . . p xi (1 − p xi ) p xi (1 − p xi ) p xi . . . ... ... . . . . . . . . . ... p xi (1 − p xi ) q − p xi (1 − p xi ) q − . . . . . . p xi , (21)where q = T − T + 1 . By combining (11), (21), we obtainthat c i ( k ) = p xi ( − ∆ s i ( k − − p xi ) c i ( k − , (22)for all k ∈ [ T + 1 , T + 1] . Meanwhile, we set c i ( k ) = 0 for all k / ∈ [ T + 1 , T + 1] , since no testing occurs. Remark 3.
The proportion of daily confirmed cases c i ( k ) in (22) consists of two terms: the first term p xi ( − ∆ s i ( k − canbe interpreted as an infected individual’s urgency in obtain-ing a test, and the second term (1 − p xi ) c i ( k − captures theunwillingness/unlikeliness of an infected individual acquiringa test. Finally, we relate the proportion of the daily number ofrecoveries, i.e., d i ( k ) , with the underlying states. In the datacollected, d i ( k ) corresponds to the change in the proportionof recovered individuals and the total number of knownactive cases A i ( k − . We assume d i ( k ) ∼ Bin (cid:16) A i ( k − N i , hγ i ( k − (cid:17) . (23)Namely, each known active case recovers with healing rate hγ i ( k − . From [30], when the number of active cases islarge, d i ( k ) is approximately equal to hγ i ( k − A i ( k − N i .The above analysis links the collected data proportionswith the underlying states of the system. If we acquire theparameter: p xi , we will be able to estimate the state systemsas follows: Definition 4.
We assume that: (cid:98) x i ( k ) = (cid:98) x i (0) , (cid:98) r i ( k ) = (cid:98) r i (0) ,and (cid:98) s i ( k ) = (cid:98) s i (0) , where (cid:98) x i (0) , (cid:98) r i (0) , (cid:98) s i (0) ∈ [0 , for all i ∈ [ n ] , k < T . Given the testing data set Ω i ( k ) collectedfrom time step T + 1 to T + 1 , according to (22) , we definethe estimated proportion of new infections at node v i as − (cid:99) ∆ s i ( k ) = c i ( k + 1) − (1 − p xi ) c i ( k ) p xi , k ∈ [ T , T ] . (24)Notice that when p xi = 1 , (24) becomes: − (cid:99) ∆ s i ( k ) = c i ( k + 1) , k ∈ [ T , T ] , which can be interpreted as: every infected individual willbe tested the day after being infected. Hence, the estimatedchange in proportion of infection on a given day k exactlyequals to the fraction of the number of positive cases on thenext day k + 1 . Moreover, we let (cid:99) ∆ r i ( k ) = 0 for k = T . Note that thefollowing equality holds from the formulation of the SIRmodel: (cid:99) ∆ s i ( k ) + (cid:99) ∆ x i ( k ) + (cid:99) ∆ r i ( k ) = 0 . We further define that (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 ) , (25) (cid:98) r i ( k ) = (cid:98) r i ( k −
1) + (cid:99) ∆ r i ( k ) , for k ∈ [ T , T ] .According to (23) and [30], the change in the proportionof recovered individuals at node v i can be inferred as (cid:99) ∆ r i ( k ) = N i d i ( k ) A i ( k − (cid:98) x i ( k − , k ∈ [ T + 1 , T ] , (26)where (cid:98) x i ( k − is calculated from (24) and (25). When A i ( k −
1) = 0 , we assume (cid:99) ∆ r i ( k ) = 0 .Therefore, if the testing data Ω i ( k ) is available over aninterval k ∈ [ T +1 , T +1] , we can estimate the states of thesystem by repetitively applying (24), (25), and (26) with theinitial conditions, i.e., (cid:98) s i (0) , (cid:98) x i (0) , and (cid:98) r i (0) , assumed forthe geometric distribution model. This addresses question (ii)in Section II. Assumption 3.
We assume that c i ( k ) = 0 for all k ∈ [ T ] ∪{ } and the initial inferred susceptible proportion is (cid:98) s i (0) . Remark 4.
When estimating the system states, we firstassume an initial condition for the system based on reality.We also assume that outside of the testing period, theproportion of positive cases collected is zero.
Proposition 1.
Under Assumption 3, the error of the infer-ence method in (24) - (26) at time k is given by | (cid:98) s i ( k ) − s i ( k ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:98) s i (0) − s i (0) − T − (cid:88) l =1 ∆ s i ( l ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (27) for all k ≥ T .Proof. From (2a), we first represent s i ( k ) by: s i ( k ) = s i (0) + k (cid:88) l =1 ∆ s i ( l ) . (28)Now, we characterize (cid:98) s i ( k ) : (cid:98) s i ( k ) = (cid:98) s i (0) + k (cid:88) l =1 (cid:99) ∆ s i ( l )= (cid:98) s i (0) + k (cid:88) l = T (cid:99) ∆ s i ( l ) (29) = (cid:98) s i (0) − c i ( k + 1) p xi − k (cid:88) l = T +1 c i ( l ) , (30)where (29) is written because − (cid:99) ∆ s i ( l ) = 0 for all l ≤ T − in (24). We acquired (30) through representing each (cid:99) ∆ s i ( l ) , l ≥ T by (24) and following Assumption 3. By applying (22), we calculate the (cid:80) kl = T +1 c i ( l ) on the R.H.S.of (30) as k (cid:88) l = T +1 c i ( l ) = − p xi k − (cid:88) l = T ∆ s i ( l ) + (1 − p xi ) k − (cid:88) l = T +1 c i ( l ) (31) = − p xi k − (cid:88) l = T ∆ s i ( l ) + (1 − p xi ) (cid:104) k (cid:88) l = T +1 c i ( l ) − c i ( k ) (cid:105) , (32)since (cid:80) k − l = T +1 c i ( l ) = (cid:80) kl = T +1 c i ( l ) − c i ( k ) . We canreorganize (32) and acquire: k (cid:88) l = T +1 c i ( l ) = − k − (cid:88) l = T ∆ s i ( l ) − − p xi p xi c i ( k ) . (33)Hence, we replace (cid:80) kl = T +1 c i ( l ) on the R.H.S. of (30)with (33) and obtain: (cid:98) s i ( k ) = (cid:98) s i (0) − c i ( k + 1) p xi + k − (cid:88) l = T ∆ s i ( l ) + 1 − p xi p xi c i ( k ) (34) = (cid:98) s i (0) + k (cid:88) l = T ∆ s i ( l ) , (35)where (35) follows from writing c i ( k + 1) in (34) as p xi ( − ∆ s i ( k )) + (1 − p xi ) c i ( k ) , using (22). Therefore, we cancalculate | (cid:98) s i ( k ) − s i ( k ) | by comparing (28) with (35) andyield the result.Prop. 1 provides an analytical expression of the estimationerror given the initial susceptible level assumed and the starttesting time. Hence, Prop. 1 solves question (iii) in Section II. Corollary 2.
In Prop. 1, if we assume that (cid:98) s i (0) = 1 , thenwe can write that (cid:98) s i ( k ) ≥ s i ( k ) , for all k ≥ T . Moreover,if we assume the inferred initial conditions correctly, i.e., (cid:98) s i (0) = s i (0) , and T = 1 , then the algorithm will estimatethe susceptible state perfectly. Remark 5.
The result of Prop. 1 consists of two parts: (cid:98) s i (0) − s i (0) and − (cid:80) T − l =1 ∆ s i ( l ) . The first componentdepends on the difference between the inferred initial sus-ceptible level and true initial susceptible level. The secondcomponent depends on the start testing date. Therefore, theaccuracy of the estimation algorithm corresponds to theestimated initial condition and how early the testing datais collected. We will explore this error via simulations in Section VI.By estimating the proportion of infected individuals ina subpopulation of a network, we are able to acquire theestimation of the infection prevalence in the whole sys-tem. These inferred states provide an understanding of theepidemic and important factors for designing eradicationschemes for infectious diseases.
V. D
ISTRIBUTED E RADICATION S TRATEGY
In this section, we propose two distributed strategies thatemploy the true states and the estimated states, respectively,and guarantee the eradication of the virus in at least expo-nential time.We propose the following healing rate to control theepidemic spread over the network: (cid:101) γ i ( k ) = s i ( k ) n (cid:88) j =1 β ij ( k ) + (cid:15) i , i ∈ [ n ] , (36)where (cid:15) i > , for each i ∈ [ n ] . This algorithm can beunderstood as boosting the healing rate of each subpopu-lation separately by providing effective medication, medicalsupplies, and/or healthcare workers. Theorem 2.
Consider the system in (2) and assume that1) ≤ h (cid:80) j β ij ( k ) < , ∀ i ∈ [ n ] and ∀ k ∈ Z ≥ ,2) B ( k ) is symmetric and irreducible ∀ k ∈ Z ≥ ,3) ∃ (cid:15) i small enough that h (cid:101) γ i ( k ) < , ∀ i ∈ [ n ] , k ∈ Z ≥ .Then the algorithm (36) guarantees GES of the set of healthystates and x ( k ) converges to with at least an exponentialrate.Proof. By substituting (36) into (2), we obtain x i ( k + 1) = x i ( k )+ h { s i ( k ) n (cid:88) j =1 β ij ( k ) x j ( k ) − [ s i ( k ) n (cid:88) j =1 β ij ( k ) + (cid:15) i ] x i ( k ) } . (37)The state transition matrix of (37) can be written as (cid:102) M ( k ) = I + h [ S ( k ) B ( k ) − ( S ( k ) diag ( B ( k ) n × )+ diag ( (cid:15) i ))] . (38)For any i, j ∈ [ n ] , j (cid:54) = i , the entries of the i -th row of (cid:102) M ( k ) are (cid:101) m ii ( k ) = 1 − h [ s i ( k ) n (cid:88) j (cid:54) = i β ij ( k ) + (cid:15) i ] , (39) (cid:101) m ij ( k ) = hs i ( k ) β i,j ( k ) , (40)which satisfies the following inequality (cid:101) m ii ( k ) + n (cid:88) j (cid:54) = i (cid:101) m ij ( k ) ≤ − h min { (cid:15) i } , ∀ i ∈ n. (41)Therefore, by Gershgorin circle theorem, the spectral radiusof (cid:102) M ( k ) is upper bounded by − h min { (cid:15) i } : ρ ( (cid:102) M ( k )) ≤ − h min { (cid:15) i } . (42)Since we have x ( k + 1) = (cid:102) M ( k ) x ( k ) and x ( k ) ≥ for all k , we can write that (cid:107) x ( k +1) (cid:107) ≤ [1 − h min { (cid:15) i } ] (cid:107) x ( k ) (cid:107) forall k . Since (cid:15) i > , ∀ i ∈ n , we obtain that, for all x i (0) ∈ [0 , n , (cid:107) x ( k ) (cid:107) ≤ [1 − h min { (cid:15) i } ] k (cid:107) x (0) (cid:107) ≤ e − kh min { (cid:15) i } (cid:107) x (0) (cid:107) , (43)where the second inequality holds by Bernoulli’s inequal-ity [38], e x = lim n →∞ (1 + xn ) n ≥ x. (44) Hence, x ( k ) converges to with an exponential rate of atleast h min { (cid:15) i } . Therefore, the set of healthy states is GES. Remark 6.
The control strategy proposed in Theorem 2can be interpreted as follows: if the healing rate of eachsubpopulation is appropriately increased according to itssusceptible proportion, for example by distributing effectivemedication, medical supplies, and/or healthcare workers toeach subpopulation, then the epidemic will be eradicatedwith at least an exponential rate. This theorem providesdecision makers insight into, given sufficient resources, howto allocate medical supplies and healthcare workers todifferent subpopulations so that the epidemic can be erad-icated quickly. Furthermore, Theorem 2 provides sufficientconditions for guaranteeing an exponentially decreasing (cid:107) x ( k ) (cid:107) for all k when the conditions apply. In other words,implementing the control strategy in Theorem 2 at full lengthwill prevent the potential upcoming waves of the epidemicin the 2-norm sense of x ( k ) . Using the estimation results from Section IV, we considerthe following healing rate: (cid:98)(cid:101) γ i ( k ) = (cid:98) s i ( k ) n (cid:88) j =1 β ij ( k ) + (cid:15) i , i ∈ [ n ] , (45)where (cid:98) s i ( k ) is the estimated susceptible rate from (25). Corollary 3.
Consider the system in (2) and assume that1) ≤ h (cid:80) j β ij ( k ) ≤ , ∀ i ∈ [ n ] and ∀ k ∈ Z ≥ ,2) B ( k ) is symmetric and irreducible ∀ k ∈ Z ≥ .3) ∃ (cid:15) i small enough that h (cid:98)(cid:101) γ i ( k ) < , (cid:98) s i (0) = 1 ∀ i ∈ [ n ] and ∀ k ∈ Z ≥ .Then the algorithm (45) guarantees GES of the set of healthystates and x ( k ) converges to with at least an exponentialrate.Proof. Similar to the proof of Theorem 2, we substitute (45)into (2) and obtain the state transition matrix for x i ( k ) : (cid:99)(cid:102) M ( k ) = I + h [ S ( k ) B ( k ) − ( (cid:98) S ( k ) diag ( B ( k ) n × )+ diag ( (cid:15) i ))] , (46)where (cid:98) S ( k ) = diag ( (cid:98) s i ( k )) . For any i, j ∈ [ n ] , j (cid:54) = i , theentries of the i -th row of (cid:99)(cid:102) M ( k ) satisfy: (cid:98)(cid:101) m ii ( k ) + n (cid:88) j (cid:54) = i (cid:98)(cid:101) m ij ( k ) ≤ − h min { (cid:15) i } , (47)since from Corollary 2 we know that when we assumethat (cid:98) s i (0) = 1 , we obtain (cid:98) s i ( k ) ≥ s i ( k ) for all i ∈ [ n ] . Consequently, by the Gershgorin circle theorem, weobtain that the spectral norm of (cid:99)(cid:102) M ( k ) is upper bounded by − h min { (cid:15) i } .Therefore, by referring to (43) and (44), we acquire thatthe set of healthy states is GES.Theorem 2 (resp. Corollary 3) has proven that given thetrue (resp. estimated) susceptible state the distributed erad-ication strategy proposed eradicates the virus with at least FIGURE
2: Graph topology in the map of the state of Indi-ana [39] analyzed and the evolution of infected proportionin each cityan exponential rate. Therefore, question (iv) from Section IIhas been addressed here.In this section, we have presented two distributed erad-ication strategies based on the true and estimated systemstates. Both strategies ensure that the SIR epidemics convergeto the sets of healthy states exponentially. We compare thetwo strategies with numerical simulations in Section VI, andstudy how a system will react if the eradication strategiesare removed too early.VI. S
IMULATIONS
In this section, we simulate a virus spreading over a staticnetwork with 5 nodes in Fig 2 to illustrate our results.The nodes are modeled after the five metropolitan areaswith a population over 150,000 in northern Indiana, U.S.:Gary (G), Lafayette (L), Indianapolis (I), Fort Wayne (F)and South Bend (S). Two nodes are neighbors if there is amajor highway connecting them. We set the initially infectedproportion to be . at node I and . at node G and 0elsewhere. The infection rates, healing rates, and the size ofeach subpopulation are static and presented in Table. I. Theevolution of the infected proportion for each city is shownin Fig. 2. β ij G L I F SG 0.08 0.15 0.24 0 0.06L 0.15 0.12 0.13 0 0I 0.24 0.13 0.25 0.05 0.04F 0 0 0.05 0.11 0.15S 0.06 0 0.04 0.14 0.09 γ i N i I: Network Parameters of Fig. 2Considering the stochastic framework, we simulate testingdata using (22) and (23), with p xi = 0 . , ∀ i ∈ { G, L, I, F, S } from T = 6 to T = 300 . The number of daily andcumulative confirmed cases and removed (recovered) casesover time at node L are shown in Fig. 3. When k ≥ ,the proportion of infected individuals at node L begins todecrease in Fig. 2, which leads to the decline of the numberof active cases in Fig. 3.We now use the method proposed in Section IV to estimatethe susceptible proportion at node I. We assume that the C L , D L k C L , D L , A L k FIGURE
3: Simulated daily and cumulative number of casesfor node L with p xi = 0 . initial condition of the recovered state is (cid:98) r I (0) = 0 . Hence,the initial infected state is written as: (cid:98) x I (0) = 1 − (cid:98) s I (0) .In Fig. 4, we plot the absolute value of the estimationerror of the susceptible state at k = 100 versus the starttesting time T and initial condition assumed, (cid:98) s I (0) . It canbe seen in Fig. 4 (top) that the estimation error increaseslinearly with the initial susceptible level assumed. Whenthe initial condition is assumed correctly for node I, witha later start testing date, the estimation error at k = 100 builds up from to r I ( k ) eventually. The increase in theestimation error with T signifies the importance of an earlytesting during an outbreak: with appropriate initial conditionsassumed, we should initiate testing as quickly as possible toimprove the accuracy of the state estimation. Meanwhile,Fig. 4 (bottom) indicates that with a later start testing date,we must assume a lower initial susceptible level accordinglyto achieve accurate estimation. The intuition behind thisfinding is that since, by Definition 4, (cid:98) s i ( k ) = (cid:98) s i (0) , forall k < T , the lower initial condition can compensate formissed tests from k ∈ [0 , T − , captured by the lastterm in (27). However, guessing (cid:98) s i (0) correctly, namely (cid:98) s i (0) = s i (0) + (cid:80) T − l =1 ∆ s i ( l ) , for T > is quite difficult.Additionally, if we assume that (cid:98) s i (0) = 1 , the estimated (cid:98) s i ( k ) is always larger than the true susceptible state in Fig. 4.The overestimation of susceptible level encourages us todesign a stronger strategy to eradicate the virus, as will beseen in the subsequent simulations.We simulate three scenarios over the network in Fig. 2with the parameters of Table. I: no control, the distributederadication strategy in (36), and the distributed eradicationstrategy utilizing estimated states in (45). The inferred stateswere produced by the algorithm in Section IV with p xi = 0 . ∀ i ∈ { G, L, I, F, S } . The average states for each scenarioare plotted in Fig. 5. It can be seen that both eradicationstrategies are able to eliminate the virus at a much higherspeed than with no control. Furthermore, when k ≥ ,the healthy states with the two eradication strategies appliedachieve higher susceptible fractions than the healthy statewithout control. We can interpret the higher susceptibleproportion as fewer individuals in the network becoming sickduring the entire outbreak. The control algorithm from (45)converges to a healthy state faster than the algorithm in (36),and both eradication strategies prevent resurgences of thevirus over the network. In Fig. 6, we remove both eradicationstrategies when k = 50 and k = 100 and do not reinstate | (cid:98) s I ( k ) − s I ( k ) | T (cid:98) s I (0) (cid:98) s I ( ) T FIGURE
4: The absolute value of susceptible state estimationerror at k = 100 with respect to start testing date and theinitial susceptible level assumed at node I, where the redpoints represent the estimation error: | (cid:98) s I ( k ) − s I ( k ) | < . .Both of the plots are illustrations of Prop. 1them. It can be seen that both of the infection curves riseup when k ≥ (resp. k ≥ ), and reach peaks beforethey slowly die out. Fig. 6 can be interpreted as removingthe allocation of resources and healthcare workers froma subpopulation too early during a pandemic, resulting inthe increase in infection level and a potential outbreak. InFig. 7, we only enforce our eradication strategies withintime interval: k ∈ [20 , and k ∈ [20 , , respectively.We can see that although the control strategies reduce theinfection level significantly, a resurgence of the outbreakoccurs immediately upon the removal of the eradicationstrategies. Hence, policy makers are suggested to enforcethe eradication strategies during the entire outbreak to avoidthe upcoming wave of epidemic.VII. C ONCLUSION
This paper studied the stability, inference and control ofdiscrete time, time-varying SIR epidemics over networks.We established the sufficient condition for GES of the setof healthy states. In addition, we proposed a stochasticframework for estimating the underlying epidemic statesfrom collected testing data. We provided analytic expressions n (cid:80) n i x i ( k ) k n (cid:80) n i s i ( k ) k FIGURE
5: Average system states over time n (cid:80) n i x i ( k ) k n (cid:80) n i x i ( k ) k FIGURE
6: Average infection proportion of the virus overtime with the eradication strategies enforced at k ∈ [0 , (left) and at k ∈ [0 , (right)for the error of the estimation algorithm. We also proposedtwo distributed control strategies that are able to eradicatethe virus in at least exponential time. The control strategiesprovide insights for decision makers on how to eliminate anongoing outbreak.In future work, we plan to study the stability andcontrol of models with more states than SIR such asSEIRS (susceptible-exposed-infected-recovered-susceptible)and SAIR (susceptible-asymptomatic-infected-recovered) asthey can possibly capture the characteristics of COVID-19 better. In our stochastic testing framework, we did notconsider the existence of inaccurate testing kits, which appearfrequently and cause confusion for policy makers. Hence,we plan to include false positive/negative test results intoour testing and estimation model and investigate the newmodel’s estimation accuracy in the future. Furthermore, weaim to apply our model on the real data to identify the systemparameters. n (cid:80) n i x i ( k ) k n (cid:80) n i x i ( k ) k FIGURE
7: Average infection proportion of the virus overtime with both of the eradication strategies imposed at k ∈ [20 , (left) and at k ∈ [20 , (right) A CKNOWLEDGMENT
The authors would like to thank Ashish Hota (IIT Kharag-pur) and Baike She (Purdue University) for useful conversa-tions related to Section IV.R
Chemical Biology & Drug Design , 2020.[3] J. A. Backer, D. Klinkenberg, and J. Wallinga, “Incubation period of2019 novel coronavirus (2019-ncov) infections among travellers fromWuhan, China, 20–28 January 2020,”
Eurosurveillance , vol. 25, no. 5,p. 2000062, 2020.[4] W.-j. Guan, Z.-y. Ni, Y. Hu, W.-h. Liang, C.-q. Ou, J.-x. He, L. Liu,H. Shan, C.-l. Lei, D. S. Hui et al. , “Clinical characteristics of coro-navirus disease 2019 in China,”
New England Journal of Medicine ,vol. 382, no. 18, pp. 1708–1720, 2020.[5] Q. Li, X. Guan, P. Wu, X. Wang, L. Zhou, Y. Tong, R. Ren, K. S.Leung, E. H. Lau, J. Y. Wong et al. , “Early Transmission Dynamicsin Wuhan, China, of Novel Coronavirus–Infected Pneumonia,”
NewEngland Journal of Medicine , 2020.[6] S. A. Lauer, K. H. Grantz, Q. Bi, F. K. Jones, Q. Zheng, H. R.Meredith, A. S. Azman, N. G. Reich, and J. Lessler, “The incubationperiod of coronavirus disease 2019 (COVID-19) from publicly re-ported confirmed cases: estimation and application,”
Annals of InternalMedicine , vol. 172, no. 9, pp. 577–582, 2020.[7] O. Byambasuren, M. Cardona, K. Bell, J. Clark, M.-L. McLaws, andP. Glasziou, “Estimating the extent of true asymptomatic COVID-19and its potential for community transmission: systematic review andmeta-analysis,”
Available at SSRN 3586675 , 2020.[8] D. Chang, G. Mo, X. Yuan, Y. Tao, X. Peng, F.-S. Wang, L. Xie,L. Sharma, C. S. Dela Cruz, and E. Qin, “Time Kinetics of ViralClearance and Resolution of Symptoms in Novel Coronavirus Infec-tion,”
American Journal of Respiratory and Critical Care Medicine ,vol. 201, no. 9, pp. 1150–1152, 2020.[9] K. Mizumoto, K. Kagaya, A. Zarebski, and G. Chowell, “Estimatingthe asymptomatic proportion of coronavirus disease 2019 (COVID-19)cases on board the diamond princess cruise ship, Yokohama, Japan,2020,”
Eurosurveillance , vol. 25, no. 10, p. 2000180, 2020.[10] A. J. Ing, C. Cocks, and J. P. Green, “COVID-19: in the Footsteps ofErnest Shackleton,”
Thorax , 2020.[11] A. Lajmanovich and J. A. Yorke, “A Deterministic Model for Gon-orrhea in a Nonhomogeneous Population,”
Mathematical Biosciences ,vol. 28, no. 3-4, pp. 221–236, 1976.[12] C. O. Dike, Z. M. Zainuddin, and I. J. Dike, “Susceptible infectedremoved epidemic model extension for efficient analysis of ebola virusdisease transmission,”
Advanced Science Letters , vol. 23, no. 9, pp.9107–9114, 2017.[13] G. C. Calafiore, C. Novara, and C. Possieri, “A modified sir modelfor the covid-19 contagion in italy,” arXiv preprint arXiv:2003.14391 ,2020.[14] K. Rock, S. Brand, J. Moir, and M. J. Keeling, “Dynamics of infectiousdiseases,”
Reports on Progress in Physics , vol. 77, no. 2, p. 026602,2014.[15] W. Mei, S. Mohagheghi, S. Zampieri, and F. Bullo, “On the dynamicsof deterministic epidemic propagation over networks,”
Annual Reviewsin Control , vol. 44, pp. 116–128, 2017.[16] R. L. Haffajee and M. M. Mello, “Thinking Globally, Acting Lo-cally—The US Response to COVID-19,”
New England Journal ofMedicine , vol. 382, no. 22, p. e75, 2020.[17] P. Van Mieghem, J. Omic, and R. Kooij, “Virus Spread in Networks,”
IEEE/ACM Trans. on Networking (TON) , vol. 17, no. 1, pp. 1–14,2009.[18] H. J. Ahn and B. Hassibi, “Global dynamics of epidemic spread overcomplex networks,” in
Proceedings of the 52nd IEEE Conference onDecision and Control , 2013, pp. 4579–4585.[19] J. Liu, P. E. Paré, E. Du, and Z. Sun, “A Networked SIS DiseaseDynamics Model with a Waterborne Pathogen,” in
Proceedings of the2019 American Control Conference (ACC) . IEEE, 2019, pp. 2735–2740. [20] N. T. Bailey et al. , The mathematical theory of infectious diseases andits applications . Charles Griffin & Company Ltd, 5a Crendon Street,High Wycombe, Bucks HP13 6LE., 1975.[21] M. Pascual and A. Dobson, “Seasonal patterns of infectious diseases,”
PLoS Medicine , vol. 2, no. 1, p. e5, 2005.[22] S. Gracy, P. Paré, H. Sandberg, and K. H. Johansson, “Analysis andDistributed Control of Periodic Epidemic Processes,”
IEEE Trans. onControl of Network Systems , 2020.[23] P. E. Paré, C. L. Beck, and A. Nedi´c, “Stability Analysis and Controlof Virus Spread over Time–Varying Networks,” in
Proceedings of the54th IEEE Conference Decision and Control , 2015, pp. 3554–3559.[24] P. E. Paré, C. L. Beck, and A. Nedi´c, “Epidemic Processes overTime-Varying Networks,”
IEEE Trans. on Control of Network Systems ,vol. 5, no. 3, pp. 1322–1334, 2018.[25] V. Bokharaie, O. Mason, and F. Wirth, “Spread of epidemics intime-dependent networks,” in
Proceedings of the 19th InternationalSymposium on Mathematical Theory of Networks and Systems–MTNS ,vol. 5, no. 9, 2010.[26] P. E. Paré, S. Gracy, H. Sandberg, and K.H. Johansson, “Data-drivendistributed mitigation strategies and analysis of mutating epidemicprocesses,” in Proceeding of the 59th IEEE Conference on Decisionand Control, arXiv preprint arXiv:2008.07317 , 2020, accepted.[27] B. A. Prakash, H. Tong, N. Valler, M. Faloutsos, and C. Faloutsos,“Virus propagation on time-varying networks: Theory and immuniza-tion algorithms,” in
Proceeding of the Joint European Conference onMachine Learning and Knowledge Discovery in Databases . Springer,2010, pp. 99–114.[28] Q. Liu, “The Threshold of a Stochastic Susceptible–Infective EpidemicModel under Regime Switching,”
Nonlinear Analysis: Hybrid Systems ,vol. 21, pp. 49–58, 2016.[29] S. Bergquist, T. Otten, and N. Sarich, “COVID-19 pandemic in theUnited States,”
Health Policy and Technology , 2020.[30] A. R. Hota, J. Godbole, P. Bhariya, and P. E. Paré, “A closed-loopframework for inference, prediction and control of SIR epidemics onnetworks,” arXiv preprint arXiv:2006.16185 , 2020.[31] E. M. Liotta, A. Batra, J. R. Clark, N. A. Shlobin, S. C. Hoffman, Z. S.Orban, and I. J. Koralnik, “Frequent neurologic manifestations andencephalopathy-associated morbidity in COVID-19 patients,”
Annalsof Clinical and Translational Neurology , vol. 7, no. 11, pp. 2221–2230,2020.[32] R. D. Mitrani, N. Dabas, and J. J. Goldberger, “COVID-19 cardiacinjury: Implications for long-term surveillance and outcomes in sur-vivors,”
Heart Rhythm , vol. 17, no. 11, pp. 1984–1990, 2020.[33] A. Rantzer, “Distributed control of positive systems,” in
Proceedingsof the 50th IEEE Conference on Decision and Control and EuropeanControl Conference , 2011, pp. 6608–6611.[34] M. Vidyasagar,
Nonlinear Systems Analysis . SIAM, 2002.[35] W. J. Rugh,
Linear System Theory . Prentice Hall Upper Saddle River,NJ, 1996, vol. 2.[36] J. An, X. Liao, T. Xiao, S. Qian, J. Yuan, H. Ye, F. Qi, C. Shen,L. Wang, Y. Liu et al. , “Clinical characteristics of recovered covid-19patients with re-detectable positive rna test,”
Annals of translationalmedicine , vol. 8, no. 17, 2020.[37] J. Abbasi, “The promise and peril of antibody testing for covid-19,”
Jama , vol. 323, no. 19, pp. 1881–1883, 2020.[38] N. L. Carothers,
Real Analysis
Matrix Analysis . Cambridge UniversityPress, 2012. A PPENDIX
A. Proof of Theorem. 1Proof.
By Assumptions 1 and 2, M ( k ) is nonnegative.Therefore, from Lemma 2, for all k ∈ Z ≥ , there exists adiagonal matrix Q ( k ) (cid:31) such that M (cid:62) ( k ) Q ( k +1) M ( k ) − Q ( k ) ≺ .Consider the following Lyapunov function V ( k, x ) = x (cid:62) Q ( k ) x . Since for all k ∈ Z ≥ , Q ( k ) is diagonal andpositive definite, it can be written that x (cid:62) Q ( k ) x > , forall x (cid:54) = . Therefore, V ( k, x ) > for all k ∈ Z ≥ , x (cid:54) = . Additionally, all the eigenvalues of Q ( k ) are real andpositive. By applying Rayleigh-Ritz Quotient Theorem [40],we obtain λ min ( Q ( k )) I ≤ Q ( k ) ≤ λ max ( Q ( k )) I, (48)which implies σ (cid:107) x (cid:107) ≤ V ( k, x ) ≤ σ (cid:107) x (cid:107) , (49)where σ = min k ∈ Z ≥ λ min ( Q ( k )) and σ =max k ∈ Z ≥ λ max ( Q ( k )) , with σ , σ > , for all k ∈ Z ≥ .Now we turn to ∆ V ( k, x ) . For x (cid:54) = 0 and for each k ∈ Z ≥ , using (3) and (6)-(7), we can write ∆ V ( k, x )= x (cid:62) ˆ M ( k ) (cid:62) Q ( k + 1) ˆ M ( k ) x − x (cid:62) Q ( k ) x = x (cid:62) [ M ( k ) (cid:62) Q ( k + 1) M ( k ) − Q ( k )] x − hx (cid:62) B (cid:62) ( k )( I − S ( k )) Q ( k + 1) M ( k ) x + h x (cid:62) B (cid:62) ( k )( I − S ( k )) Q ( k + 1)( I − S ( k )) B ( k ) x. (50)Note that the second and third term of (50) can be reorga-nized as x (cid:62) [ − hB (cid:62) ( k )( I − S ( k )) Q ( k + 1) M ( k )+ h B (cid:62) ( k )( I − S ( k )) Q ( k + 1)( I − S ( k )) B ( k )] x = x (cid:62) { hB (cid:62) ( k )( I − S ( k )) Q ( k + 1)[ − M ( k ) + h ( I − S ( k )) B ( k )] } x = x (cid:62) { hB (cid:62) ( k )( I − S ( k )) Q ( k + 1)[ − I − h Γ( k )) − h ( I + S ( k )) B ( k )] } x ≤ , (51)where the last equality follows from (6), and the inequalityfollows from Assumptions 1 and 2 and Lemma 1. Thus, byapplying (51) into (50), we obtain that ∆ V ( k, x ) ≤ x (cid:62) [ M ( k ) (cid:62) Q ( k + 1) M ( k ) − Q ( k )] x. (52)From Lemma 2, we know that [ M ( k ) (cid:62) Q ( k +1) M ( k ) − Q ( k )] is negative definite, and [ Q ( k ) − M ( k ) (cid:62) Q ( k + 1) M ( k )] ispositive definite. Therefore, we obtain λ max [ M ( k ) (cid:62) Q ( k +1) M ( k ) − Q ( k )] = − λ min [ Q ( k ) − M ( k ) (cid:62) Q ( k + 1) M ( k )] .By applying Rayleigh-Ritz Quotient Theorem we can write ∆ V ( k, x ) ≤ − σ (cid:107) x (cid:107) , (53)where σ = max k ∈ Z ≥ λ min [ Q ( k ) − M ( k ) (cid:62) Q ( k + 1) M ( k )] ,with σ > for all k ∈ Z ≥ .Therefore, from (49), (53) and Lemma 3 the set of healthystates of (2) is GES, proving the theorem. B. Proof of Corollary 1Proof.
From Lemma 4, (49), and (53), the rate of conver-gence is upper bounded by (cid:113) − σ σ . The next step is toshow that the rate is well defined, which is (cid:113) − σ σ ∈ [0 , .Since σ > and σ > , we only need to prove that σ ≥ σ . Note that for all k ∈ Z ≥ , Q ( k ) and M ( k ) (cid:62) Q ( k +1) M ( k ) are both symmetric. Therefore, by applying Weyl’s inequali-ties from [40] to [ Q ( k ) − M ( k ) (cid:62) Q ( k + 1) M ( k )] , we obtainthat, for all k ∈ Z ≥ , λ i [ Q ( k ) − M ( k ) (cid:62) Q ( k + 1) M ( k )] ≤ λ i ( Q ( k )) − λ i [ M ( k ) (cid:62) Q ( k + 1) M ( k )] . (54)We compare the LHS of (54) with σ and the RHS of (54)with σ to yield σ ≤ λ i [ Q ( k ) − M ( k ) (cid:62) Q ( k + 1) M ( k )] , (55) σ ≥ λ i ( Q ( k )) − λ i [ M ( k ) (cid:62) Q ( k + 1) M ( k )] , (56)where (56) holds because M ( k ) (cid:62) Q ( k + 1) M ( k ) is positivesemidefinite for all k ∈ Z ≥ . Hence, σ ≥ σ3