key: cord-0540870-9t373sxm
authors: Wang, Dan; Liu, Ji; Par'e, Philip E.; Chen, Wei; Qiu, Li; Beck, Carolyn L.; Bacsar, Tamer
title: Controlling a Networked SIS Model via a Single Input over Undirected Graphs
date: 2020-04-23
journal: nan
DOI: nan
sha: a13801c5ae0bb67e57222d2f428db72a046ce2be
doc_id: 540870
cord_uid: 9t373sxm

This paper formulates and studies the problem of controlling a networked SIS model using a single input in which the network structure is described by a connected undirected graph. A necessary and sufficient condition on the values of curing and infection rates for the healthy state to be exponentially stable is obtained via the analysis of signed Laplacians when the control input is the curing budget of a single agent. In the case when the healthy state is stabilizable, an explicit expression for the minimum curing budget is provided. The utility of the algorithm is demonstrated using a simulation over a network of cities in the northeastern United States.

Mathematical models of virus spread have been studied for centuries (Bernoulli, 1760) . These models have offered interesting insight into how spread processes appear in real world systems (Paré et al., 2020) . This paper focuses on SIS (susceptible-infected-susceptible) models that were first introduced by Kermack and McKendrick (1932) . The idea is that every member of the population is either healthy (susceptible) or infected, and transitions occur from being infected to being healthy according to some curing rate and from being healthy to infected according to some infection rate.

The original SIS model (Kermack and McKendrick, 1932) had assumed that the population is completely connected and has instantaneous mixing, modeling the population as two groups, a susceptible group and an infected group.

Researchers have extended the model to capture more realistic structures, including nontrivial graphs in the mathematical models (Lajmanovich and Yorke, 1976; Fall et al., 2007; Mieghem et al., 2009; Khanafer et al., 2016) . There are two interpretations of these models and their states: 1) that each node is a subpopulation and the state is a proportion of infected individuals in the corresponding subpopulation (Lajmanovich and Yorke, 1976; Fall et al., 2007), or 2) that each node is a single individual and the state is the probability of that individual being infected (Mieghem et al., 2009) . Various extensions of these models have been proposed to include time-varying graph structures (Paré et al., 2015 (Paré et al., , 2018 , and competing viruses on layered networks (Santos et al., 2015; Watkins et al., 2018; Liu et al., 2019) .

There have been various control techniques proposed for networked SIS models; for example, Torres et al. (2017) ; Wan et al. (2007) ; Khanafer et al. (2016) ; Preciado et al. (2014) ; Watkins et al. (2018) ; Mai et al. (2018) ; Drakopoulos et al. (2014) , to name a few. Drakopoulos et al. (2014) proposed a dynamic curing policy for containment of a contagion process modeled by a continuous Markov process. It was shown that if the CutWidth of the underlying graph is sublinear in the number of nodes, the proposed policy achieved sublinear expected time to extinction. Borgs et al. (2010) studied the epidemic threshold when the antidote was distributed non-uniformly. Wan et al. (2007 Wan et al. ( , 2008 proposed and implemented distributed control techniques for setting healing rates and quarantine protocols on a severe acute respiratory syndrome (SARS) simulation model. The work of Khanafer et al. (2016) has proposed an antidote control technique. Preciado et al. (2014) and Watkins et al. (2018) developed an optimal vaccination control technique using geometric programming ideas. Mai et al. (2018) use distributed optimization to solve several formulations of curing resource allocation problems. The proposed algorithms by Mai et al. (2018) , while efficient and mostly distributed, all appeal to an additional, independent consensus process to calculate a piece of centralized information, and thus require a synchronous stopping time across the network. None of the other proposed algorithms in the literature are distributed. In fact, it has been shown by Liu et al. (2019) that a certain type of distributed feedback controller can never stabilize the healthy state.

With the preceding discussions in mind, we are interested in the following minimum control budget problem. Suppose that the networked SIS system is unstable at the healthy state, which implies that the network will ultimately converge to an epidemic state as long as at least one agent gets infected (Lajmanovich and Yorke, 1976; Fall et al., 2007; Mieghem et al., 2009) , and that we can select a fixed set of agents to boost their curing rates. This technique gives a fixed, additive control input to their curing rates that is in place for all time. Then, is it possible to stabilize the healthy state? If possible, what is the minimum total control budget needed to stabilize the healthy state? These problems turn out to be challenging and, to the best of our knowledge, there is no existing algorithm or analytic solution.

In this paper, we study a special case of the problems identified above, in which only one control input is allowed. This special case is a good starting point for several reasons. Consider the interpretation of the model where each node is a subpopulation; if a government has enough resources to treat only one subpopulation, this approach can show if the disease can be eradicated given that only one subpopulation can be treated, and if so, which is the best subpopulation to treat, and how much effort will be required to eradicate the disease. An analogous interpretation can also be made for a computer network to be protected by equipping virus patches on only one computer. In the case when the underlying graph is undirected, we obtain a necessary and sufficient condition for the healthy state to be exponentially stable via the mathematical tool of signed Laplacians. When the condition is satisfied, i.e., when the healthy state is stabilizable, we provide an explicit expression for the minimum control budget.

Notation: For any positive integer n, we use [n] to denote the set {1, 2, . . . , n}. We use 0 and 1 to denote the vectors whose entries all equal 0 and 1, respectively, and I to denote the identity matrix, while the dimensions of the vectors and matrices are to be understood from the context. For any two real vectors a, b ∈ IR n , we write

. We use M † to denote the Moore-Penrose inverse of a matrix M . The eigenvalues of a symmetric matrix A ∈ R n×n are denoted by λ 1 (A) ≤ · · · ≤ λ n (A). We write A 0 and A 0 when the matrix A is (symmetric) positive definite and positive semidefinite, respectively.

The continuous-time single-virus SIS model over an nagent network (after mean-field approximation) is as follows (Lajmanovich and Yorke, 1976) :ẋ

where x i (t) denotes the probability of agent i to be infected, δ i is the nonnegative curing rate of agent i, and β ij is the nonnegative infection rate between agents i and j, whose value equals zero if they are not neighbors. The neighbor relationships among the n agents are described by an undirected n-vertex graph G whose vertices represent the agents and edges depict the neighbor relationships. We assume that β ij = β ji for all i, j ∈ [n].

The n equations in (1) can be combined into one equation in a state-space form aṡ

where x(t) is the state vector in IR n whose ith entry is x i (t), ∆ is the n × n diagonal matrix whose ith diagonal entry is δ i , B is the n × n matrix whose ijth entry is β ij , and X(t) is the n × n diagonal matrix whose ith diagonal entry is x i (t). Since we have assumed that β ij = β ji for all i, j ∈ [n], B is symmetric. It is also assumed that B is an irreducible matrix, which implies that the underlying graph G is connected.

It is easy to see that 0 is an equilibrium of system (2). We call it the healthy state. It has been shown in Lajmanovich and Yorke (1976) ; Mieghem et al. (2009); Fall et al. (2007) that the healthy state is the unique equilibrium of the system if ∆ − B 0, and if ∆ − B 0 does not hold, the healthy state is unstable, and there will be a unique epidemic equilibrium x * 0, which is almost globally stable. More can be said. The following result is a variant of Proposition III.1 in Nowzari et al. (2017) . Proposition 1. The healthy state x = 0 of system (2) is exponentially stable if, and only if, ∆ − B 0.

With this in mind, we formulate the following control problem which aims to stabilize the healthy state using a single control input.

Suppose that the necessary and sufficient condition in Proposition 1 is not met, and we have a control budget to increase the curing rate of a single agent i, i.e,

We seek to answer the questions: what is the necessary and sufficient condition to stabilize the healthy state? When the healthy state is stabilizable with such a single control input, what is the minimum control budget? Formally stated we have: minimize

where U i is the n × n diagonal matrix whose ith diagonal entry equals u i and the other diagonal entries all equal zero. In the following, we will present machinery that enables solving this problem in an optimal and efficient manner.

Mathematically, the problems in (3) and (4) amount to making a symmetric matrix positive definite by changing one of its diagonal entries with minimal effort. We are aware that this problem, in its more general form of changing multiple diagonal entries, has been explored in the literature based on linear matrix inequalities (LMI) (Bhatia, 2009 ), which provides numerical solutions but not analytical ones.

The recent works (Torres et al., 2017; Torres and Roy, 2018) studied the problem of minimizing the dominant eigenvalue of ∆ + U i − B when a subset of nodes can be controlled, subject to a given total control budget, with u i ≥ 0 in the first paper and u i ∈ R in the second one.

Algorithms have been proposed for optimal subset design. This problem is different from the one we consider in this paper. Here, we wish to minimize the total control budget such that ∆ + U i − B 0. Moreover, we want to obtain an analytic solution for the minimum control budget when the healthy state is stabilizable. We achieve these goals for the special case of changing a single diagonal entry by exploiting the structure of signed graphs.

Let D i be the n×n diagonal matrix whose diagonal entries are

Let V i be an n × n diagonal matrix whose ith diagonal entry equals (δ i − n j=1 β ij ) and all the other diagonal entries equal zero. Define

Note that 1 D i 1 = n j=1 d j . By construction, the row sums of L i all equal zero.

The following theorem provides an explicit lower bound on the control input u i that stabilizes the healthy state. Theorem 1. Suppose that u i is the single control input of agent i. The healthy state of system (3) is exponentially stable if, and only if, L i 0 with a simple eigenvalue at zero and

where e i is the (n + 1)-dimensional unit vector with the ith entry being one and others zero.

It is worth emphasizing that L i is independent of u i , as is the right side of (6). The matrix L i 0 and 0 is a simple eigenvalue with eigenspace span{1}. Thus L † i 0 and 0 is a simple eigenvalue with the eigenspace span{1}, which will be shown in the next section. Since (e i − e n+1 ) / ∈ span{1}, the last term in (6) is positive. If the value of the right side of (6) is negative, it implies that the healthy state is exponentially stable without any control input, i.e., ∆ − B 0. Given any i, by Theorem 1, since L i is independent of u i , we can first check whether the curing resources on the other agents are sufficient to stabilize the healthy state or not, i.e., whether the condition that L i 0 with a simple eigenvalue at zero holds or not. Subsequently, if that condition is met, then we can design a control input, appealing to (6), to make the healthy state exponentially stable. Corollary 1. Suppose that u is the single control input on any one agent. The healthy state of system (3) is stabilizable if, and only if, there exists at least one i ∈ [n] such that L i 0 with a simple eigenvalue at zero. Let J be the set of such indices in [n] . If the condition holds, i.e., J = ∅, the minimum control input to stabilize the healthy state is given by

where u i satisfies (6).

The corollary is a direct consequence of Theorem 1 and provides the machinery to explicitly solve the problem posed in (4).

To prove the main results stated in the previous section, we need the following preliminaries on results for Laplacian matrices and their variants.

Consider a weighted undirected graph G = (V, E) with no self-loops, which consists of a set of nodes V = [n] and a set of edges E. We associate with each edge (i, j) ∈ E a real-valued weight w ij . The associated Laplacian matrix L = [l ij ] n×n is defined by

It is clear that L is a symmetric matrix whose row sums all equal zero, and thus zero is an eigenvalue of L. It is worth emphasizing that we allow negative weights w ij . In the case where all weights w ij are nonnegative, it is well known that L is positive semidefinite. Moreover, zero is a simple eigenvalue with eigenspace span{1} if and only if G is connected (Fiedler, 1973) . In the case when there exists at least one negative weight w ij , L is sometimes called a signed Laplacian, whose property of positive semidefiniteness has been recently investigated in Zelazo and Bürger (2014) , Chen et al. (2016b) , and Chen et al. (2016a) . With this work in mind, any symmetric matrix whose row sums all equal zero can be regarded as a Laplacian matrix whose corresponding weighted undirected graph may contain negative weights.

There is a useful factorization of a Laplacian matrix L. Suppose that the corresponding graph G has m edges. Let W = diag{w 1 , w 2 , . . . , w m } denote the m × m diagonal matrix where the diagonal entries are the edge weights, i.e., for each edge ε k ∈ E, w k = w ij for (i, j) = ε k . The order of the m edges can be arbitrary. In addition, we assign an arbitrary orientation to each edge, i.e., for each edge ε k ∈ E, we set an arbitrary endpoint as the head and the other one as the tail. The oriented incidence matrix

otherwise. Then, the Laplacian matrix L can be factorized as L = BW B .

(8) It is worth noting that while the incidence matrix B depends on the choice of orientations, the Laplacian matrix L does not. In the case when both positive and negative weights exist, we define G + and G − as the spanning subgraphs of G which consist of all positive-and negativeweighted edges, and we call them positive and negative subgraphs, respectively. In this case, L can also be written in the following form:

where B + and B − are incidence matrices corresponding to the positive and negative subgraphs, respectively, and W + and W − are the diagonal matrices whose diagonal entries are the weights of the edges in the positive and negative subgraphs, respectively. From (8), it is easy to see that B + W + B + and B − W − B − are the Laplacian matrices of G + and G − , respectively. In fact, it can be verified that if G can be partitioned into p spanning subgraphs, G 1 , G 2 , . . . , G p , where all p subgraphs have the same set of vertices, if their edge sets are mutually disjoint, and if the union of their edge sets equals the edge set of G, then the Laplacian L can be written as

where B i is the incidence matrix of G i and W i is the diagonal matrix whose diagonal entries are the weights of the edges in G i . It can be seen that each B i W i B i equals the Laplacian matrix of G i .

We will later make use of the following lemmas. Lemma 1. A Laplacian matrix L satisfies L 0 with a simple eigenvalue at zero and associated eigenspace span{1} if, and only if, L † 0 with a simple eigenvalue at zero and associated eigenspace span{1}.

Proof: Consider the Schur decomposition of L as

where U 1/ √ n1 is a unitary matrix and S is a diagonal matrix containing all the nonzero eigenvalues of L. By (Campbell and Meyer, 2009 , Theorem 1.2.1), we have

which is the Schur decomposition of L † . Therefore, the statement holds. Lemma 2. Let L be an n × n Laplacian matrix and Q be a matrix in R n×(n−1) such that Q Q = I and Q 1 = 0.

Then, L 0 with simple zero eigenvalue if, and only if, Q LQ 0.

Proof: Since Q Q = I and Q 1 = 0, it follows that

is a unitary matrix. Note that

Then, the spectrum of L consists of zero and the eigenvalues of Q LQ, from which the lemma follows. Lemma 3. Let L 0 be an n × n Laplacian matrix with a simple zero eigenvalue and Q be a matrix in R n×(n−1) such that Q Q = I and Q 1 = 0. Then, L † = Q(Q LQ) −1 Q .

Proof: It can be seen that Q 1/ √ n1 is a unitary matrix. By Theorem 1.2.1 in Campbell and Meyer (2009) and Lemma 2, we have

= Q(Q LQ) −1 Q , which completes the proof.

Consider a weighted undirected graph G = (V, E) which consists of a set of nodes V = [n] and a set of edges including self-loops. Associate with each edge (i, j) ∈ E a real-valued weight w ij . The weight of the self-loop of node i is denoted by w ii ∈ R. The associated loopy Laplacian H (Dorfler and Bullo, 2013) is defined as H = L + D, where L is the Laplacian matrix as in (7), and D = diag{w 11 , . . . , w nn } is a diagonal matrix whose diagonal entries are the self-loop weights. When D = 0, i.e., at least one of w ii is nonzero, H is called a strictly loopy Laplacian (Dorfler and Bullo, 2013) . Clearly, H and L are symmetric matrices.

Given a strictly loopy Laplacian, we introduce a grounded node with index n + 1, and connect node i ∈ [n] to the grounded node when w ii = 0; the self-loops are then eliminated. This process gives us a loopless grapĥ G (i.e., there are no self-loops in the graph) with either positive or negative weights, and its signed Laplacian L ∈ R (n+1)×(n+1) . The augmented LaplacianL is defined asL

which is also a symmetric matrix.

Lemma 4. Consider a strictly loopy Laplacian matrix H ∈ R n×n and its corresponding augmented Laplacian matrix L. Then, H 0 if, and only if,L 0 with a simple zero eigenvalue.

Proof: We first prove the sufficiency. Suppose thatL 0 and zero is a simple eigenvalue. SinceL1 = 0 and zero is a simple eigenvalue, the eigenspace of the zero eigenvalue is span{1}. Let x ∈ R n be any nonzero vector and set y = [x 0] . Since y / ∈ span{1}, from (11), we have x Hx = y L y > 0, which implies that H is positive definite.

We now prove the necessity. Suppose that H 0. Denote the eigenvalues of H andL by λ 1 (H) ≤ · · · ≤ λ n (H) and λ 1 (L) ≤ · · · ≤ λ n+1 (L), respectively. Since H is a principal submatrix ofL, by the Cauchy Interlace Theorem (see Theorem 4.3.8 in Horn and Johnson (1990) ), we have λ 1 (L) ≤ λ 1 (H) ≤ λ 2 (L) ≤ · · · ≤ λ n (H) ≤ λ n+1 (L). Thus, λ i (L) > 0 for all i ∈ {2, . . . , n + 1}. Since the row sums ofL are all equal to zero,L has a zero eigenvalue, which implies that λ 1 (L) = 0. This completes the proof.

To prove the theorem, we will need the following lemma. To proceed, we rewrite the matrix ∆ + U i − B as follows. LetD i be the diagonal matrix whose diagonal entries arê

Then, ∆+U i −B −D i is a matrix whose row sums all equal zero, and thus can be regarded as a Laplacian matrix.

can be treated as a strictly loopy Laplacian of a weighted graph in which the weight on each edge (i, j) is β ij and the weight of the self-loop of each node j isd j . To proceed, let

Note that 1 D i 1 = n j=1d j . It is easy to check that all row sums ofL i are equal to zero. Note,L i is the augmented Laplacian for ∆ + U i − B. Specifically, we have introduced a grounded node with index n + 1 to the existing graph G, as explained in Section 4.2, to create a loopless grapĥ G, and thus removing the self-loops from the graph. The Fig. 1 . Plot of the northeast train routes: we ignore the routes that lead out of this subset of cities.

corresponding signed LaplacianL i ∈ R (n+1)×(n+1) , given in (14), has possibly both positive and negative weights.

Note that the difference betweenD i and D i only lies in the ith diagonal entry. It is not hard to verify that L i given in (5) is the Laplacian matrix of the graph obtained by removing the edge (i, n + 1) fromĜ.

From Proposition 1 and Lemma 4, it is sufficient to show thatL i 0 with a simple zero eigenvalue if and only if L i 0 with a simple zero eigenvalue and (6) holds.

Let Q be a matrix in R (n+1)×n such that

is a unitary matrix. Consequently, Q Q = I and Q 1 = 0. Then, from Lemma 2,L i 0 with a simple zero eigenvalue if and only if Q L i Q 0. SinceL i is a Laplacian matrix whose corresponding graphĜ's edges are all positiveweighted except for the edge (i, n + 1) whose sign depends on the value of u i . From (10),

in which (e i − e n+1 )d i (e i − e n+1 ) equals the Laplacian of the spanning subgraph ofĜ whose edge set has only one edge (i, n + 1). Thus,L i 0 with a simple zero eigenvalue if and only if

Ifd i > 0, thenĜ is a connected graph whose edge weights are all positive, which implies thatL i 0 with a simple eigenvalue at zero. Ifd i = 0, thenL i = L i . From Lemma 1, it can be seen that L i 0 is equivalent to (6), which implies that the statementL i 0 with a simple zero eigenvalue is equivalent to the statement L i 0 with a simple zero eigenvalue and (6) holds. We therefore consider the case whend i < 0.

From item 1) in Lemma 5, (15) is equivalent to

From item 2) in Lemma 5, (16) is equivalent to the following two conditions:

The first condition (17) is equivalent to that L i 0 with a simple zero eigenvalue because of Lemma 2. The second condition (18) is equivalent to (6) due to Lemma 3. This then completes the proof. Remark 1. A key element in the above proof is showing that the condition thatL i is positive semidefinite with a simple zero eigenvalue is equivalent to the condition that L i is positive semidefinite with a simple zero eigenvalue and (6) holds. Alternative proofs of such equivalence are available by adapting the techniques in the proofs of Theorem III.3 in Zelazo and Bürger (2014) and Theorem 3 in Song et al. (2018) . Still, the new proof provided as above is much more concise compared to the existing techniques in the literature.

For the simulations we consider a virus spreading among the northeast of the United States. The graph structure is comprised of nine cities, or subpopulations, in the northeast (see Figure 1 ). We make the simplifying assumption that there is a virus which spreads among subpopulations, grouped by cities, according to the passenger train network (ignoring spread among air and car travel). The connectivity is determined by the Amtrak Rail Service between these cities. See Figure 1 for a plot of the train routes. The corresponding infection rate matrix is 

We set the healing rate to δ i = 3.5 for all i. The state of the system is the proportion of infected individuals in the corresponding city. Red indicates completely infected (x i = 1) and blue indicates completely healthy (x i = 0).

The necessary and sufficient condition from Proposition 1 is not satisfied; the minimum eigenvalue of ∆ − B is −0.3134. Therefore, with no control effort the system converges to the epidemic equilibrium

x * = [.063 .034 .034 .026 .104 .080 .031 .090 .070] , meaning between 3-10% of the populations of the cities are infected with the virus at steady state. Now the question becomes the following: if we are limited by our infrastructure and can treat only one city, which one is the best one and how much treatment is required? Implementing the algorithm from Corollary 1 to solve Problem (4), answers these questions. It turns out that the only node that meets the base requirement of L i 0 with a simple eigenvalue at zero is i = 5, that is, New York City (node 5). The minimum amount of control input that eradicates the virus is u 5 = 1.30. Including this control input in the system as in (3) results in the system converging to the healthy state, depicted in Figure 2 .

From this example, we discovered, as would be expected since there is only one control input, that the virus cannot be too strong for a system to be stabilizable using this t x i (t) Fig. 2 . The healthy equilibrium of the northeast system (left, blue indicates healthy) and trajectories (right) with infection rates and connectivity given in (19), the healing rate is δ i = 3.5 for all i, and the control input u 5 = 1.30.

approach. If δ i is much smaller than 3.5 in the system considered here, then there is not an i such that L i 0.

We have proposed a control technique that employs one control input on a single node to mitigate the spread of a virus in a network, that can drive the system to the healthy state. A necessary and sufficient condition on the values of curing and infection rates for the healthy state to be exponentially stable has been obtained. In the case when the healthy state is stabilizable, an explicit expression for the minimum curing budget has been provided. The derivation of the controller used a loopy Laplacian approach, employing a new, more concise proof technique than what can be found in the literature. We have shown the utility of the proposed controller via simulation.

For future work we would like to explore, using a similar control technique, situations with more than one control input. We would also like to extend the technique to be distributed across the network. Finally, we would like to implement these techniques in some real spread applications, to eradicate virus from a real system.

Essai dune nouvelle analyse de la mortalité causée par la petite vérole et des avantages de linoculation pour la prévenir

Positive Definite Matrices

How to distribute antidote to control epidemics. Random Structures &

Generalized Inverses of Linear Transformations. SIAM

Characterizing the positive semidefiniteness of signed Laplacians via effective resistances

On the definiteness of graph Laplacians with negative weights: geometrical and passivity-based approaches

Kron reduction of graphs with applications to electrical networks

An efficient curing policy for epidemics on graphs

A Course in Robust Control Theory: A Convex Approach

Epidemiological models and Lyapunov functions

Algebraic connectivity of graphs

Matrix Analysis

Contributions to the mathematical theory of epidemics. II. The problem of endemicity

Stability of epidemic models over directed graphs: a positive systems approach

A deterministic model for gonorrhea in a nonhomogeneous population

Analysis and control of a continuoustime bi-virus model

Distributed algorithm for suppressing epidemic spread in networks

Virus spread in networks

Optimal resource allocation for control of networked epidemic models

Stability analysis and control of virus spread over time-varying networks

Epidemic processes over time-varying networks

Analysis, estimation, and validation of discretetime epidemic processes

Optimal resource allocation for network protection against spreading processes

Bi-virus SIS epidemics over networks: qualitative analysis

Network-based analysis of small-disturbance angle stability of power systems

Dominant eigenvalue minimization with trace preserving diagonal perturbation: Subset design problem

Sparse resource allocation for linear network spread dynamics

Network design problems for controlling virus spread

Designing spatially heterogeneous strategies for control of virus spread

Optimal resource allocation for competitive spreading processes on bilayer networks

On the definiteness of the weighted Laplacian and its connection to effective resistance