key: cord-0559758-kcm9pqws
authors: Zhang, Ciyuan; Leung, Humphrey; Butler, Brooks; Par'e, Philip. E.
title: Estimation and Distributed Eradication of SIR Epidemics on Networks
date: 2021-02-24
journal: nan
DOI: nan
sha: dea22b6b5e9cdc621e8be9ef208d1241a3e569ed
doc_id: 559758
cord_uid: kcm9pqws

This work examines the discrete-time networked SIR (susceptible-infected-recovered) epidemic model, where the infection and recovery parameters may be time-varying. We provide a sufficient condition for the SIR model to converge to the set of healthy states exponentially. We propose a stochastic framework to estimate the system states from observed testing data and provide an analytic expression for the error of the estimation algorithm. Employing the estimated and the true system states, we provide two novel eradication strategies that guarantee at least exponential convergence to the set of healthy states. We illustrate the results via simulations over northern Indiana, USA.

As of February 2021, the COVID-19 virus has claimed 2.4 million lives and infected 110 million individuals worldwide [1] . Lack of effective treatments, high contagion rates [2] , long incubation periods [3] - [6] , and asymptomatic cases [7] - [10] pose significant challenges in containing and eradicating pandemics. Recent pandemics, including gonorrhea [11] , Ebola [12] , and COVID-19 [13] , have accelerated the development of infection models. The main goal of epidemic model development is to identify conditions to eradicate the pathogen, and leverage the knowledge of these conditions to design mitigation strategies. Various infection models have been proposed, based on characteristics of individual pathogens, and studied in the literature, including susceptible-infected-susceptible (SIS), susceptibleinfected-removed (SIR), and susceptible-infected-removedsusceptible (SIRS) [14] , [15] . In this paper, we focus on the SIR epidemic model. We aim to expand on the SIR model, by exploring mutating viruses over networks, estimation of the underlying states, and distributed eradication strategies.

The patchwork response to COVID-19 [16] gives rise to susceptible community subpopulations, with heterogeneous time-varying factors not previously explored by the SIR model. Extensions on the SIS model, studied in [17] - [19] , augment the compartmental epidemic models originated in [20] to include interactions between subpopulations of susceptible communities. Additionally, various advanced epidemic models consider time-varying factors [21] - [28] . In this paper, we establish sufficient conditions for the set of healthy states of a networked time-varying SIR model to be globally exponentially stable (GES). These equilibrium states are not unique, as the final susceptible and removed states *Ciyuan Zhang, Humphrey Leung, Brooks Butler, and Philip. E. Paré are with the School of Electrical and Computer Engineering at Purdue University. Emails: {zhan3375, leung61, brooksbutler, philpare}@purdue.edu. *This work was funded in part by the C3.ai Digital Transformation Institute sponsored by C3.ai Inc. and the Microsoft Corporation and in part by the National Science Foundation, grants NSF-CNS #2028738 and NSF-ECCS #2032258.

are dependent on the initial conditions and the time-varying infectious and healing parameters.

The delay in onset of COVID-19 symptoms [3] - [6] , large asymptomatic populations estimated between 17 − 81% [7]- [10] , and delay in test results [29] compromise the ability for accurate estimation of current infection states. An estimation algorithm that incorporated a constant delay between the change in infection proportion and testing data was introduced in [30] . They studied the inference problem by using a Bayesian approach. Inspired by the delay characterization suggested in [30] , we propose a stochastic delay to model the unpredictability of the COVID-19 virus and testing strategies. We have developed methods for estimating the underlying epidemic states from testing data with a delay sampled from a geometric distribution, which cannot be completely filtered by the method suggested in [30] . The geometric delay model accounts for the stochastic effect of individuals failing to get tested immediately after exposure. We study the aggregated effect of each individual delay on the trajectory of confirmed cases and devise a method for estimating the underlying epidemic states of an SIR model from these delayed measurements. We also investigate the proposed method's estimation error, which provides insights for achieving an accurate estimation of the system states. We then employ this more realistic estimation to strategically eradicate a disease.

As proven in [30] , the SIR epidemic model converges to a healthy state, however, an exponential convergence is not shown. Combining the modeling and inference approach allows us to develop two distributed control strategies capable of eradicating epidemic spread exponentially, at an equilibrium with a higher proportion of susceptible population. Decreasing the removed (recovered) and increasing the susceptibility proportion is of exceptional importance for the COVID-19 pandemic, as long term and severe health complications have been documented in the recovered populations, including impaired cognition [31] and damage to cardiac tissue [32] . Our main result shows that by applying the estimated and the true susceptible states of each node, our proposed eradication strategies will guarantee global exponential stability of a healthy state of the overall network.

We summarize the main contributions of this paper as follows:

• We establish sufficient conditions for global exponential stability of the set of healthy states; see Theorem 1. • We propose a stochastic framework which estimates the trajectories of the system states of the networked SIR model from testing data. Furthermore, we provide analytical expressions for the error of the estimation algorithm we propose; see Prop. 1. • We propose two distributed eradication strategies for adjusting healing rates, one that is based on the true system states and the other one is based on the inferred system states. Both methods guarantee that the virus is eradicated within exponential time; see Theorem 2 and Corollary 3.

We organize this paper as follows: Section II lays down some basic assumptions and restates the well-known SIR model in the networked fashion, and it presents the main problems studied in this paper. Section III first recalls preliminary results that are essential for stability analysis and then discusses the sufficient conditions for global exponential stability of a healthy state of the networked time-varying SIR models. Section IV covers the proposed techniques of estimating hidden epidemics states with the stochastic delay of tested individuals and testing data. Section V covers the two distributed control strategies which ensure that the system converges to a healthy state in at least exponential time. Section VI illustrates the results from Section IV and V with numerical simulations. Finally, in Section VII, we summarize the main conclusions of this paper and discuss potential future directions.

We denote the set of real numbers, the non-negative integers, and the positive integers as R, Z ≥0 , and Z ≥1 , respectively. For any positive integer n, we have [n] = {1, 2, ..., n}. The spectral radius of a matrix A ∈ R n×n is ρ(A). A diagonal matrix is denoted as diag(·). The transpose of a vector x ∈ R n is x . The Euclidean norm is denoted by · . We use I to denote identity matrix. We use 0 and 1 to denote the vectors whose entries all equal 0 and 1, respectively. The dimensions of the vectors are determined by context. Given a matrix A, A 0 (resp. A 0) indicates that A is positive definite (resp. positive semidefinite), whereas A ≺ 0 (resp. A 0) indicates that A is negative definite (resp. negative semidefinite). Let G = (V, E) denote a graph or network where V = {v 1 , v 2 , ..., v n } is the set of subpopulations, and E ⊆ V × V is the set of edges. We denote the expectation of a random variable as IE[·].

Consider a time-varying epidemic network of n subpopulations, where the size of subpopulation v i is N i ∈ Z >0 , and the infection rates and healing rates could be timevarying. 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 rate of node v i at time t. The proportions of the subpopulation at node v i which are susceptible, infected, and recovered at time t are denoted by s i (t), x i (t) and r i (t), respectively. The deterministic continuous-time evolution of the SIR epidemic is given bẏ

We now state the discrete-time SIR epidemic dynamics obtained through Euler discretization of (1). For a small sampling time h > 0, the discrete-time evolution of the SIR epidemic is given by

Equation (2b) can be rewritten as

where S(k) = diag(s(k)), B(k) is the matrix with (i, j)th entry β ij (k), and Γ(k) = diag(γ i (k)). The spread of a virus over a network can be captured using a graph G = (V, E),

is the set of directed edges. We make the following assumptions in order for the system in (2) to be well defined. We have the following result which shares the same idea as the time-invariant model, proven in [30] . 

Definition 1. We define the set of healthy states of (2) as

. Given a network that is infected by a virus, our goal is to guarantee that each subpopulation v i converges to the set of healthy states in exponential time regardless of the initial conditions of the each subpopulations. We now officially state the questions being studied in this paper:

(i) For the system with dynamics given in (3), under what condition is the set of healthy states, i.e., x(k) = 0, global exponentially stable (GES)?

(ii) Given the testing data, how can the stochastic framework be constructed to estimate the susceptible, infected, and recovered proportions, denoted by s i (k), x i (k), and r i (k), respectively, for each subpopulation v i in the network? (iii) What is the estimation error of the stochastic framework that we proposed? (iv) Given the knowledge of the conditions that ensure GES of a healthy state, i.e., x(k) = 0 and s i (k) inferred from testing data, how can we devise dynamic control algorithms which apply new healing rates γ i (k) to each agent in (3) so that the epidemic is eradicated with a faster rate of convergence than the rate of exponential?

This section presents conditions that ensure global exponential stability of the set of healthy states. First, we introduce some preliminaries and then we present our main analysis results.

In this subsection, we recall results that are crucial for understanding the rest of the paper.

Suppose that M is a nonnegative matrix which satisfies ρ(M ) < 1. Then there exists a diagonal matrix P 0 such that M P M − P ≺ 0.

Consider a system described as follows:

Definition 2. An equilibrium point of (4) is GES is there exist positive constants α and ω, with 0 ≤ ω < 1, such that

We recall a sufficient condition for GES of an equilibrium of (4) from [34] .

, and ∀x(k 0 ) ∈ R n , then x(k) = 0 is a globally exponential stable equilibrium of (4). Note that a healthy state of the system in (2) is GES if Assumptions 1 and 2 hold and the condition in Definition 2 and Lemma 3 are satisfied for all x(k 0 ) ∈ [0, 1] n , since this is the domain where the model is well defined.

In this subsection, we present sufficient conditions for the global exponential stability of the set of healthy states of the system. We find the conditions by analyzing the spectral radius of the state transition matrix of (2b). We define

Notice thatM (k) is the state transition matrix of (2b) and it can be written that

We use M (k) and (8) to illustrate the sufficient conditions for the GES of the set of healthy states in the subsequent theorem.

Theorem 1. Given Assumptions 1 and 2, suppose for all

then the set of healthy states of (2) is GES.

Proof. See Appendix.

Corollary 1. Under the assumptions of Theorem 1, the rate of convergence to a healthy state is upper bounded by an exponential rate of

Proof. See Appendix.

Remark 1. Notice that in Theorem 1, sup k∈Z ≥0 ρ(M (k)) < 1 is the key condition which ensures the set of healthy states of (2) is GES. We can interprete ρ(M ) in the context of epidemiology as the basic reproduction number of the virus over the network. Particularly, Theorem 1 affirms that given the time-varying parameters of the network satisfy the condition provided, the mutating virus will exponentially converge to the set of healthy states.

In this section, we found the conditions that ensure exponential convergence to the set of healthy states of (2b) in Theorem 1. This answers question (i) in Section II. The stability condition can help the policymakers reallocate the medical resources, staff etc. which leads to modifying the parameters in (2) so that the spreading of the virus stops completely. One of the other factors that will assist in decision making is the COVID-19 observed testing data.

In this section, we study how to estimate the epidemic states (s(k), x(k), r(k)) from testing data in order to design a feedback controller in the following section. One of the challenges of estimating the underlying system states is that the testing data on a given day does not capture the new infections on the same day. Instead, the testing data is a delayed representation of the change in the system. Characterizing the delay of each individual is difficult, because the delay is determined by numerous factors such as the incubation period of COVID-19, the duration of obtaining test results, the willingness of each individual to get tested, etc. Therefore, we propose a stochastic framework in this section to capture the factors which cause the testing delay. In our discrete-time model, we assume that τ i ∈ Z ≥0 . We model the testing delay of each infected individual τ i as two aggregate components to represent the uncertainty in the testing process:

where η i ∈ Z ≥0 is a constant and Y i is sampled from a discrete time random variable whose measurable space is Z ≥0 .

Remark 2. In (10), the constant component η i can be interpreted as the length of time needed to acquire testing results. The random variable can be interpreted as the incubation period and/or the amount of time that it takes an individual to get tested after becoming infected.

First, we denote the set of estimated system states for

is the number of confirmed cases at time k, and D i (k) represents the number of removed (recovered) cases at time k. In addition, the cumulative number of confirmed and removed cases at node v i are written as C i (k) = k j=0 C i (j) and D i (k) = k j=0 D i (j), respectively. Therefore, the number of active cases is calculated by

Ni as the proportion of daily confirmed cases and removal, respectively. Note that

The estimation procedure is illustrated in Fig. 1 .

We then study how to relate c i (k) to the underlying states. We define a vector space Π T1 as the space of all the proportions of daily number of confirmed cases from time step k = T 1 to time step k = T 2 + 1. We define Ξ T1 as the vector of all the decreases in the proportion of susceptible individuals, −∆s i (k), from time step k = T 1 to time step k = T 2 + 1. We denote Φ(T 1 , T 2 ) as the transfer matrix which results in

where

which depends on the SIR dynamics, the testing strategies, and the delay. When the testing delay is a constant, i.e.,

When the delay η i = 0, the transfer matrix Φ(T 1 ,

We now propose a stochastic testing framework to capture the delay between when an individual is infected and when they receive a positive test result. We first let η i = 0 in (10), without the loss of generality. Furthermore, we assume that each infected individual at node v i has an equal probability p x i ∈ (0, 1] of receiving a diagnostic test each day starting from the day after they are infected. Therefore, we model Y i in (10) as a random variable following the geometric distribution, with the probability of an infected individual acquiring a positive test δ days after infection being:

for δ ∈ Z ≥1 . The geometric distribution of the testing delay models the number of days before an infected individual obtains a diagnostic test which represents the incubation period of COVID-19 and/or the unwillingness of each individual getting a test. We assume that the delay of each infected individual's distribution is i.i.d. (independent and identically distributed) from others. Furthermore, we assume that an infected individual can be tested only once. Based on [36] and [37] , we assume that even if an individual recovers from COVID-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 binary random variable X i (ν) with X i (ν) = 1 (resp. X i (ν) = 0) if a randomly chosen individual from subpopulation v i became infected at time ν. It can be written that

where, from (2)

We define the binary random variable T i (µ, δ), with T i (µ, δ) = 1 if a randomly chosen individual acquired a positive test at time µ and was infected δ days before µ. Now we rewrite ν, in (15) , as µ − δ. From (14) , the conditional probability P (T i (µ, δ) = 1|X i (µ − δ) = 1) is given by the geometric probability mass function (pmf) p x i (1−p x i ) δ−1 and represents the probability of an infected individual acquiring a positive test specifically δ days after being infected. Hence, the joint pmf of the two random variables X i (µ−δ), T i (µ, δ) is written as:

which is interpreted as the probability that a randomly chosen individual became infected at time µ − δ and acquired a positive test at time µ, where µ − δ, µ ∈ [T 1 , T 2 ]. Therefore, the joint pmf P Xi,Ti (µ − δ, µ) is calculated as:

since we assume that the test results are accurate. Similarly,

Let W i (µ) be the marginal distribution of T i (µ, δ) over the set of feasible delays, δ, with its pmf being the probability of a random individual acquiring a positive test at time µ:

by combining (17) and (18) . Therefore, the number of confirmed cases at time k is calculated as

where (19) (20) for all k ∈ [T 1 + 1, T 2 + 1], the transfer matrix Φ(T 1 , T 2 ) in (11) is written as

where q = T 2 − T 1 + 1. By combining (11), (21), we obtain that

for all k ∈ [T 1 + 1, T 2 + 1]. Meanwhile, we set c i (k) = 0 for all k / ∈ [T 1 + 1, T 2 + 1], since no testing occurs. Finally, we relate the proportion of the daily number of recoveries, i.e., d i (k), with the underlying states. In the data collected, d i (k) corresponds to the change in the proportion of recovered individuals and the total number of known active cases A i (k − 1). We assume

Namely, each known active case recovers with healing rate hγ i (k − 1). From [30] , when the number of active cases is large, d i (k) is approximately equal to hγi(k−1)Ai(k−1) Ni . The above analysis links the collected data proportions with the underlying states of the system. If we acquire the parameter: p x i , we will be able to estimate the state systems as follows:

Definition 4. We assume that:

, k < T 1 . Given the testing data set Ω i (k) collected from time step T 1 + 1 to T 2 + 1, according to (22) , we define the estimated proportion of new infections at node v i as

Notice that when p x i = 1, (24) becomes: − ∆s i (k) = c i (k + 1), k ∈ [T 1 , T 2 ], which can be interpreted as: every infected individual will be tested the day after being infected. Hence, the estimated change in proportion of infection on a given day k exactly equals to the fraction of the number of positive cases on the next day k + 1.

Moreover, we let ∆r i (k) = 0 for k = T 1 . Note that the following equality holds from the formulation of the SIR model:

We further define that

According to (23) and [30] , the change in the proportion of recovered individuals at node v i can be inferred as

where x i (k − 1) is calculated from (24) and (25) . When A i (k − 1) = 0, we assume ∆r i (k) = 0. Therefore, if the testing data Ω i (k) is available over an interval k ∈ [T 1 +1, T 2 +1], we can estimate the states of the system by repetitively applying (24), (25) , and (26) with the initial conditions, i.e., s i (0), x i (0), and r i (0), assumed for the geometric distribution model. This addresses question (ii) in Section II.

Assumption 3. We assume that c i (k) = 0 for all k ∈ [T 1 ] ∪ {0} and the initial inferred susceptible proportion is s i (0). Remark 4. When estimating the system states, we first assume an initial condition for the system based on reality. We also assume that outside of the testing period, the proportion of positive cases collected is zero. 

for all k ≥ T 1 .

Proof. From (2a), we first represent s i (k) by:

Now, we characterize s i (k):

where (29) 

We can reorganize (32) and acquire:

Hence, we replace k l=T1+1 c i (l) on the R.H.S. of (30) with (33) and obtain:

where (35) follows from writing c i (k + 1) in (34) as p x i (−∆s i (k)) + (1 − p x i )c i (k), using (22) . Therefore, we can calculate | s i (k) − s i (k)| by comparing (28) with (35) and yield the result.

Prop. 1 provides an analytical expression of the estimation error given the initial susceptible level assumed and the start testing time. Hence, Prop. 1 solves question (iii) in Section II. Corollary 2. In Prop. 1, if we assume that s i (0) = 1, then we can write that s i (k) ≥ s i (k), for all k ≥ T 1 . Moreover, if we assume the inferred initial conditions correctly, i.e., s i (0) = s i (0), and T 1 = 1, then the algorithm will estimate the susceptible state perfectly. ∆s i (l). The first component depends on the difference between the inferred initial susceptible level and true initial susceptible level. The second component depends on the start testing date. Therefore, the accuracy of the estimation algorithm corresponds to the estimated initial condition and how early the testing data is collected.

We will explore this error via simulations in Section VI.

By estimating the proportion of infected individuals in a subpopulation of a network, we are able to acquire the estimation of the infection prevalence in the whole system. These inferred states provide an understanding of the epidemic and important factors for designing eradication schemes for infectious diseases.

In this section, we propose two distributed strategies that employ the true states and the estimated states, respectively, and guarantee the eradication of the virus in at least exponential time.

We propose the following healing rate to control the epidemic spread over the network:

where i > 0, for each i ∈ [n]. This algorithm can be understood as boosting the healing rate of each subpopulation separately by providing effective medication, medical supplies, and/or healthcare workers.

Theorem 2. Consider the system in (2) and assume that 1) 0 ≤ h j β ij (k) < 1, ∀i ∈ [n] and ∀k ∈ Z ≥0 , 2) B(k) is symmetric and irreducible ∀k ∈ Z ≥0 , 3) ∃ i small enough that h γ i (k) < 1, ∀i ∈ [n], k ∈ Z ≥0 . Then the algorithm (36) guarantees GES of the set of healthy states and x(k) converges to 0 with at least an exponential rate.

Proof. By substituting (36) into (2), we obtain

The state transition matrix of (37) can be written as

(38) For any i, j ∈ [n], j = i, the entries of the i-th row of M (k) are

which satisfies the following inequality

Therefore, by Gershgorin circle theorem, the spectral radius of M (k) is upper bounded by 1 − h min{ i }:

Since we have x(k + 1) = M (k)x(k) and x(k) ≥ 0 for all k, we can write that x(k +1) ≤ [1−h min{ i }] x(k) for all k. Since i > 0, ∀i ∈ n, we obtain that, for all x i (0) ∈ [0, 1] n ,

where the second inequality holds by Bernoulli's inequality [38] ,

Hence, x(k) converges to 0 with an exponential rate of at least h min{ i }. Therefore, the set of healthy states is GES.

Remark 6. The control strategy proposed in Theorem 2 can be interpreted as follows: if the healing rate of each subpopulation is appropriately increased according to its susceptible proportion, for example by distributing effective medication, medical supplies, and/or healthcare workers to each subpopulation, then the epidemic will be eradicated with at least an exponential rate. This theorem provides decision makers insight into, given sufficient resources, how to allocate medical supplies and healthcare workers to different subpopulations so that the epidemic can be eradicated quickly. Furthermore, Theorem 2 provides sufficient conditions for guaranteeing an exponentially decreasing x(k) for all k when the conditions apply. In other words, implementing the control strategy in Theorem 2 at full length will prevent the potential upcoming waves of the epidemic in the 2-norm sense of x(k).

Using the estimation results from Section IV, we consider the following healing rate:

where s i (k) is the estimated susceptible rate from (25) .

Corollary 3. Consider the system in (2) and assume that 1) 0 ≤ h j β ij (k) ≤ 1 , ∀i ∈ [n] and ∀k ∈ Z ≥0 , 2) B(k) is symmetric and irreducible ∀k ∈ Z ≥0 .

3) ∃ i small enough that h γ i (k) < 1, s i (0) = 1 ∀i ∈ [n] and ∀k ∈ Z ≥0 . Then the algorithm (45) guarantees GES of the set of healthy states and x(k) converges to 0 with at least an exponential rate.

Proof. Similar to the proof of Theorem 2, we substitute (45) into (2) and obtain the state transition matrix for x i (k):

where S(k) = diag( s i (k)). For any i, j ∈ [n], j = i, the entries of the i-th row of M (k) satisfy:

since from Corollary 2 we know that when we assume that s i (0) = 1, we obtain s i (k) ≥ s i (k) for all i ∈ [n]. Consequently, by the Gershgorin circle theorem, we obtain that the spectral norm of M (k) is upper bounded by 1 − h min{ i }. Therefore, by referring to (43) and (44), we acquire that the set of healthy states is GES.

Theorem 2 (resp. Corollary 3) has proven that given the true (resp. estimated) susceptible state the distributed eradication strategy proposed eradicates the virus with at least an exponential rate. Therefore, question (iv) from Section II has been addressed here.

In this section, we have presented two distributed eradication strategies based on the true and estimated system states. Both strategies ensure that the SIR epidemics converge to the sets of healthy states exponentially. We compare the two strategies with numerical simulations in Section VI, and study how a system will react if the eradication strategies are removed too early.

In this section, we simulate a virus spreading over a static network with 5 nodes in Fig 2 to illustrate our results. The nodes are modeled after the five metropolitan areas with 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 a major highway connecting them. We set the initially infected proportion to be 0.02 at node I and 0.01 at node G and 0 elsewhere. The infection rates, healing rates, and the size of each subpopulation are static and presented in Table. I. The evolution of the infected proportion for each city is shown in Fig. 2 Fig. 2 Considering the stochastic framework, we simulate testing data using (22) and (23) , with p x i = 0.2, ∀i ∈ {G, L, I, F, S} from T 1 = 6 to T 2 = 300. The number of daily and cumulative confirmed cases and removed (recovered) cases over time at node L are shown in Fig. 3 . When k ≥ 80, the proportion of infected individuals at node L begins to decrease in Fig. 2 , which leads to the decline of the number of active cases in Fig. 3 .

We now use the method proposed in Section IV to estimate the susceptible proportion at node I. We assume that the In Fig. 4 , we plot the absolute value of the estimation error of the susceptible state at k = 100 versus the start testing time T 1 and initial condition assumed, s I (0). It can be seen in Fig. 4 (top) that the estimation error increases linearly with the initial susceptible level assumed. When the initial condition is assumed correctly for node I, with a later start testing date, the estimation error at k = 100 builds up from 0 to r I (k) eventually. The increase in the estimation error with T 1 signifies the importance of an early testing during an outbreak: with appropriate initial conditions assumed, we should initiate testing as quickly as possible to improve 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 accordingly to achieve accurate estimation. The intuition behind this finding is that since, by Definition 4, s i (k) = s i (0), for all k < T 1 , the lower initial condition can compensate for missed tests from k ∈ [0, T 1 − 1], captured by the last term in (27) . However, guessing s i (0) correctly, namely s i (0) = s i (0) + T1−1 l=1 ∆s i (l), for T 1 > 0 is quite difficult. Additionally, if we assume that s i (0) = 1, the estimated s i (k) is always larger than the true susceptible state in Fig. 4 . The overestimation of susceptible level encourages us to design a stronger strategy to eradicate the virus, as will be seen in the subsequent simulations.

We simulate three scenarios over the network in Fig. 2 with the parameters of Table. I: no control, the distributed eradication strategy in (36) , and the distributed eradication strategy utilizing estimated states in (45). The inferred states were produced by the algorithm in Section IV with p x i = 0.5 ∀i ∈ {G, L, I, F, S}. The average states for each scenario are plotted in Fig. 5 . It can be seen that both eradication strategies are able to eliminate the virus at a much higher speed than with no control. Furthermore, when k ≥ 200, the healthy states with the two eradication strategies applied achieve higher susceptible fractions than the healthy state without control. We can interpret the higher susceptible proportion as fewer individuals in the network becoming sick during 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 the virus over the network. In Fig. 6 , we remove both eradication strategies when k = 50 and k = 100 and do not reinstate them. It can be seen that both of the infection curves rise up when k ≥ 50 (resp. k ≥ 100), and reach peaks before they slowly die out. Fig. 6 can be interpreted as removing the allocation of resources and healthcare workers from a subpopulation too early during a pandemic, resulting in the increase in infection level and a potential outbreak. In Fig. 7 , we only enforce our eradication strategies within time interval: k ∈ [20, 50] and k ∈ [20, 150], respectively. We can see that although the control strategies reduce the infection level significantly, a resurgence of the outbreak occurs immediately upon the removal of the eradication strategies. Hence, policy makers are suggested to enforce the eradication strategies during the entire outbreak to avoid the upcoming wave of epidemic.

This paper studied the stability, inference and control of discrete time, time-varying SIR epidemics over networks. We established the sufficient condition for GES of the set of healthy states. In addition, we proposed a stochastic framework for estimating the underlying epidemic states from collected testing data. We provided analytic expressions for the error of the estimation algorithm. We also proposed two distributed control strategies that are able to eradicate the virus in at least exponential time. The control strategies provide insights for decision makers on how to eliminate an ongoing outbreak.

In future work, we plan to study the stability and control of models with more states than SIR such as SEIRS (susceptible-exposed-infected-recovered-susceptible) and SAIR (susceptible-asymptomatic-infected-recovered) as they can possibly capture the characteristics of COVID-19 better. In our stochastic testing framework, we did not consider the existence of inaccurate testing kits, which appear frequently and cause confusion for policy makers. Hence, we plan to include false positive/negative test results into our testing and estimation model and investigate the new model's estimation accuracy in the future. Furthermore, we aim to apply our model on the real data to identify the system parameters. x = 0. Additionally, all the eigenvalues of Q(k) are real and positive. By applying Rayleigh-Ritz Quotient Theorem [40] , we obtain λ min (Q(k))I ≤ Q(k) ≤ λ max (Q(k))I,

which implies

where σ 1 = min k∈Z ≥0 λ min (Q(k)) and σ 2 = max k∈Z ≥0 λ max (Q(k)), with σ 1 , σ 2 > 0, for all k ∈ Z ≥0 . Now we turn to ∆V (k, x). For x = 0 and for each k ∈ Z ≥0 , using (3) and (6) 

Note that the second and third term of (50) can be reorganized as 

where σ 3 = max k∈Z ≥0 λ min [Q(k) − M (k) Q(k + 1)M (k)], with σ 3 > 0 for all k ∈ Z ≥0 . Therefore, from (49), (53) and Lemma 3 the set of healthy states of (2) is GES, proving the theorem.

Proof. From Lemma 4, (49), and (53), the rate of convergence is upper bounded by 1 − σ3 σ2 . The next step is to show that the rate is well defined, which is 1 − σ3 σ2 ∈ [0, 1). Since σ 2 > 0 and σ 3 > 0, we only need to prove that σ 2 ≥ σ 3 .

Note that for all k ∈ Z ≥0 , Q(k) and M (k) Q(k+1)M (k) are both symmetric. Therefore, by applying Weyl's inequalities from [40] to [Q(k) − M (k) Q(k + 1)M (k)], we obtain that, for all k ∈ Z ≥0 , (54)

We compare the LHS of (54) with σ 3 and the RHS of (54) with σ 2 to yield

where (56) holds because M (k) Q(k + 1)M (k) is positive semidefinite for all k ∈ Z ≥0 . Hence, σ 2 ≥ σ 3 and the rate of convergence is well defined.

World Health Organization (WHO)

The recent challenges of highly contagious COVID-19, causing respiratory infections: Symptoms, diagnosis, transmission, possible vaccines, animal models, and immunotherapy

Incubation period of 2019 novel coronavirus (2019-ncov) infections among travellers from Wuhan, China

Clinical characteristics of coronavirus disease 2019 in China

Early Transmission Dynamics in Wuhan, China, of Novel Coronavirus-Infected Pneumonia

The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: estimation and application

Estimating the extent of true asymptomatic COVID-19 and its potential for community transmission: systematic review and meta-analysis

Time Kinetics of Viral Clearance and Resolution of Symptoms in Novel Coronavirus Infection

Estimating the asymptomatic proportion of coronavirus disease 2019 (COVID-19) cases on board the diamond princess cruise ship

COVID-19: in the Footsteps of Ernest Shackleton

A Deterministic Model for Gonorrhea in a Nonhomogeneous Population

Susceptible infected removed epidemic model extension for efficient analysis of ebola virus disease transmission

A modified sir model for the covid-19 contagion in italy

Dynamics of infectious diseases

On the dynamics of deterministic epidemic propagation over networks

Thinking Globally, Acting Locally-The US Response to COVID-19

Virus Spread in Networks

Global dynamics of epidemic spread over complex networks

A Networked SIS Disease Dynamics Model with a Waterborne Pathogen

The mathematical theory of infectious diseases and its applications. Charles Griffin & Company Ltd, 5a Crendon Street, High Wycombe, Bucks HP13 6LE

Seasonal patterns of infectious diseases

Analysis and Distributed Control of Periodic Epidemic Processes

Stability Analysis and Control of Virus Spread over Time-Varying Networks

Epidemic Processes over Time-Varying Networks

Spread of epidemics in time-dependent networks

Data-driven distributed mitigation strategies and analysis of mutating epidemic processes

Virus propagation on time-varying networks: Theory and immunization algorithms

The Threshold of a Stochastic Susceptible-Infective Epidemic Model under Regime Switching

COVID-19 pandemic in the United States

A closed-loop framework for inference, prediction and control of SIR epidemics on networks

Frequent neurologic manifestations and encephalopathy-associated morbidity in COVID-19 patients

COVID-19 cardiac injury: Implications for long-term surveillance and outcomes in survivors

Distributed control of positive systems

Nonlinear Systems Analysis. SIAM

Linear System Theory

Clinical characteristics of recovered covid-19 patients with re-detectable positive rna test

The promise and peril of antibody testing for covid-19

39] Indiana state of USA

Matrix Analysis

Lemma 2, for all k ∈ Z ≥0 , there exists a diagonal matrix Q(k) 0 such that M (k)Q(k +1)M (k)− Q(k) ≺ 0. Consider the following Lyapunov function V (k, x) = x Q(k)x. Since for all k ∈ Z ≥0 , Q(k) is diagonal and positive definite

The authors would like to thank Ashish Hota (IIT Kharagpur) and Baike She (Purdue University) for useful conversations related to Section IV.