Optimal curing policy for epidemic spreading over a community network with heterogeneous population

# Optimal curing policy for epidemic spreading over a community network with heterogeneous population

Stefania Ottaviano
Fondazione Bruno Kessler, Via alla Cascata 56/d, 38123 Povo (Trento), Italy
Mathematics Department, University of Trento, Via Sommarive 14, 38123 Povo (Trento), Italy
Francesco De Pellegrini
Fondazione Bruno Kessler, Via alla Cascata 56/d, 38123 Povo (Trento), Italy
Stefano Bonaccorsi
Mathematics Department, University of Trento, Via Sommarive 14, 38123 Povo (Trento), Italy
and
Piet Van Mieghem
EEMCS, Delft University of Technology, Mekelweg 4 2628 CD Delft, Netherlands
###### Abstract

The design of an efficient curing policy, able to stem an epidemic process at an affordable cost, has to account for the structure of the population contact network supporting the contagious process. Thus, we tackle the problem of allocating recovery resources among the population, at the lowest cost possible to prevent the epidemic from persisting indefinitely in the network. Specifically, we analyze a susceptible–infected–susceptible epidemic process spreading over a weighted graph, by means of a first-order mean-field approximation. First, we describe the influence of the contact network on the dynamics of the epidemics among a heterogeneous population, that is possibly divided into communities. For the case of a community network, our investigation relies on the graph-theoretical notion of equitable partition; we show that the epidemic threshold, a key measure of the network robustness against epidemic spreading, can be determined using a lower-dimensional dynamical system. Exploiting the computation of the epidemic threshold, we determine a cost-optimal curing policy by solving a convex minimization problem, which possesses a reduced dimension in the case of a community network. Lastly, we consider a two-level optimal curing problem, for which an algorithm is designed with a polynomial time complexity in the network size. heterogeneous SIS model, community network, graph spectra, equitable partitions, convex optimization.
2000 Math Subject Classification: 82C20, 34A34, 34D23, 90C22, 90C25

The diffusion and persistence of infectious diseases depend on complex interactions between individual units (namely people, cities, countries, etc.), the characteristics of a disease and, possibly, on the applied control policies. The last ones are aimed at arresting disease transmission or render the infection prevalence as low as possible.

Epidemic models have been used to describe a wide range of other phenomena as well, like social behaviors, diffusion of information, computer viruses, etc., indeed, although the basic mechanisms of these phenomena can be different, often their dynamical behavior can be described by the same type of equations [1]. One of the main objectives, in all these domains, is to gain insight into how the spreading process transmits and to identify the most effective strategies in order to prevent and control them.

In controlling the diffusion processes, the structure of the contact network plays a crucial role. In particular, several contact networks appear organized into communities. In this framework, a uniform control strategy not always represents the most effective way to reduce the infection rate, the number of affected individuals or the time of extinction [2, 3, 4]. Furthermore, curing costs may vary from node to node. In the case of community networks, curing costs may vary depending on the features of the specific community where curing controls are applied.

Thus, by taking into account the topology of a community network, in this work we want to determine a cost-optimal distribution of resources, that is able to prevent the disease from persisting indefinitely in the population. The non-uniform distribution of resources aims to control, in a cost-optimal way, the level of the nodes local curing rates. Increasing the curing rates of, e.g., some selected communities, is reflected into speeding up their detection capabilities and treatments (or, into installing better virus scan software, in the case of computer viruses) [2].

The problem of designing strategies to stop spreading processes in networks has been largely tackled. Though, in this context, to the best of our knowledge, very few works have described how to exploit the community structures in order to formulate an optimization problem for resources allocation, with lower complexity. Based on the theory of contact processes, Borgs et al. [4] characterize the optimal distribution of a fixed amount of antidote in a given network. Gourdin et al. [5] and Sahnen et al. [6] take advantage of the -intertwined approximation [7] to analyze and control the spread of a SIS epidemic model. The same mean-field approach is adopted by Preciado et al in [8], where the authors propose a semidefinite programming (SDP) approach for optimal network immunization. Cost-optimal vaccine allocation in arbitrary undirected networks are obtained via the minimization of a vaccination cost function which depends on infection rates. In [9], the same authors specialize some specific instances of optimal network protection problem, via Geometric Programming techniques, to weighted, directed, strongly (and not necessarily strongly) connected networks to compute the cost and speed optimal allocation of vaccines and/or antidotes. In Sec. Optimal curing policy for epidemic spreading over a community network with heterogeneous population, we analyze more in detail the differences between their and our approach. Enyioha et al. [10] propose a distributed solution to the vaccine and antidote allocation problem to contain an epidemic outbreak in the absence of a central social planner. Each node locally computes its optimum investment in vaccine and antidotes needed to globally contain the spread of an outbreak, via local exchange of information with neighbors. Drakopoulos et. al [11, 12] consider the propagation of an epidemic process over a network and study the problem of dynamically allocating a fixed curing budget to the nodes of the graph. The objective is to minimize the expected extinction time of the epidemics. In the case of bounded degree graphs, they provide a lower bound on the expected time to extinction under any such dynamic allocation policy.

Compared to previous works on optimal curing policy, we are interested, particularly, in leveraging the subdivision of the population into communities. The motivation comes from the fact that community structures are a relevant non-trivial topological feature of complex networks. Community structures have been identified as a typical feature of social networks, tightly connected groups of nodes in the World Wide Web usually correspond to pages on common topics, whereas in the biology framework, e.g., in cellular and genetic networks, communities may relate to functional modules [13]. Consequently, in many practical situations, it appears reasonable to consider curing policies which apply (i.e., per hospital, school, village, or city, etc,…), rather than policies which apply per individual unit.

In particular, we consider a continuous-time susceptible–infected–susceptible (SIS) epidemic process, where an individual can be repeatedly infected, recover and yet be infected again. An input weighted graph captures the interaction between individuals and communities, where the heterogeneity, and possible asymmetries in the contagiousness, are caught by edge-dependent weights.

Our investigation, on a population divided into communities, has been based on the graph-theoretical notion of            equi [14, 15, 16]. A network with an equitable partition of its node set posses some interesting symmetry properties; we will use the word “symmetry” to refer to a certain structural regularity of the graph connectivity [16]. We take advantage of the notion of equitable partitions for providing curing policies, diversified for communities, capable to lead the system to extinction, at the minimum cost. In this context, our main goal is that we are able to formulate a convex cost minimization problem with a reduced dimension, with respect to the general case, where curing policies are providing for each node.

Spatial inhomogeneity has been incorporated in other models to study the epidemic control [17, 18, 19], however not much effort has been made to explore inhomogeneous control strategies within this kind of models [2]. The problem of an inhomogeneous allocation of limited resources for a multi-group model, has been studied, instead, in [2] by Wan et al. The aim of the authors is to maximize the speed at which the virus is eliminated. Thus, considering a discrete-time epidemic model, they tackle the problem to minimize the dominant eigenvalue of a system matrix, subject to limiting constraints on some system parameters to be controlled. Compared to their formulation, we want to allocate resources to the communities, sufficient to lead the system towards the epidemic extinction, with the aim of minimizing a certain cost function. Moreover, in the cited paper, individuals transmit the disease through homogeneous mixing within their own group, as well as interactions with individuals in other groups, like in the usual metapopulation models. Such model is only a specific case of our network model, in fact, by the notion of equitable partition, we go beyond the full mesh assumption, within the communities, as well as outside, thus providing results for wider possible scenario (see Section Optimal curing policy for epidemic spreading over a community network with heterogeneous population).

The paper is organized as follows. In Sec. Optimal curing policy for epidemic spreading over a community network with heterogeneous population, first, we review some background concepts for epidemic processes on networks in the homogeneous setting and the related mean-field approximation adopted in the paper. Then, we provide the adaptation of the model to the heterogeneous setting and we report on the analysis of the global dynamics of the epidemic process. This allows to recognize the stability modulus of a matrix, encoding for the network structure and for the parameters of the model, as the critical value separating an absorbing phase, from an endemic phase. Leveraging on this result, in Sec. Optimal curing policy for epidemic spreading over a community network with heterogeneous population, we present the cost-optimal resource allocation problem. We use a semidefinite approach to formulate our optimization problem for the case of arbitrary undirected graphs with symmetric weights. Then, we show that this approach can be extended to some kind of not symmetric weighted networks, those whose adjacency matrix is diagonally symmetrizable. Moreover, for the case of a general directed weighted graph, we provide a suboptimal solution. In Sec Optimal curing policy for epidemic spreading over a community network with heterogeneous population we consider the case when a contact network is shaped by an existing community network. First, we extend the results in [20] (related to equitable partitions in the case of undirected graphs), considering equitable partitions for directed weighted graphs. Specifically, we show how a certain kind of structure regularity, in a directed weighted graph, influences the system of differential equations that solve for the evolution in time of the approximated infection probability of each node. Then, we exploit such regularities in graph connectivity for reducing the original dynamical system to a lower dimensional one. By supposing that different curing rates can be chosen depending on the community network structure and that they can be optimized for a certain cost function, the latter system is used to reduce the dimension of our optimization problem. In the last part of the work, we propose a two-level optimal curing problem, that is, we have a two-dimensional curing policy, suitable, e.g., when the population can be divided in two categories (young and elders, male and female, etc,…). This kind of situation fits well certain networks with equitable partition, such as, e.g., bipartite graphs and interconnected star networks. In this case we provide a scalable bisection algorithm, that yields an -approximation of the optimal solution, in polynomial time in the input size. Finally, we carry out some numerical experiments. Proofs not included in previous sections for better readability are placed in Section Optimal curing policy for epidemic spreading over a community network with heterogeneous population.

In this section, we report some background concepts and new tools that we will use later to find a cost-optimal curing policy.

Let us consider a SIS epidemic process spreading over a simple undirected graph , with edge set and node (vertex) set . The order of , denoted by , is the cardinality of . The edge set of consists of unordered pairs , with , and . Connectivity of graph is conveniently expressed by the symmetric adjacency matrix .

The viral state of a node , at time , is described by a Bernoulli random variable , where we set , if is healthy and , if is infected. Every node at time is either infected with probability or healthy (but susceptible) with probability . In the homogeneous setting, the recovery process is a Poisson process with rate , and the infection process is a per link Poisson process where the infection rate between an healthy and an infected node is . All the infection and recovery processes are independent. The SIS process, developing a graph with nodes, can be modeled as a continuous-time Markov process with states, covering all possible combinations in which nodes can be infected [7]. The probability of the process of being in a certain state can be uniquely determined by the Kolmogorov’s differential equations (i.e. a system of linear differential equations). However, the number of equations increases exponentially with the number of nodes; this poses several limitations in order to determines the set of solutions even for small network order. Hence, often, it is necessary to derive models that are an approximation of the exact original one [7, 21].

In this work, we consider a first order mean-field approximation (NIMFA), proposed by Van Mieghem et. al. in [7]. Basically, NIMFA replaces the original linear differential equations by non-linear differential equations representing the time-change of the infection probability of each node.

Epidemic threshold. For a network with finite order , the exact SIS Markov process will always converge towards its unique absorbing state, that is the zero-state where all nodes are healthy. However, the process shows a phase transition behavior: indeed, there is a critical value of the effective spreading rate , whereby if is lesser than , the initial infection dies out quickly. Conversely, for larger than , the infection spreading can last very long in any sufficiently large network [22, 23, 24]. The regime of persistent infection ( ), called the metastable or quasi-stationary state, is reached rapidly given an initial set of infected nodes and can persist for very long time [24]. In support of this, numerical simulations of SIS processes reveal that, even for fairly small networks and when , the overall-healthy state is only reached after an unrealistically long time. Hence, the indication of the model is that, in the case of real networks, one should expect that above the epidemic threshold the extinction of epidemics is hardly ever attained [22, 23]. For this reason, the literature is mainly concerned with establishing the value of the epidemic threshold, being a key measure of the robustness against epidemic spreading.
In the homogeneous setting, NIMFA determines the epidemic threshold for the effective spreading rate as , where is the spectral radius of the adjacency matrix , (see [7, 25]). When the only equilibrium of the NIMFA system is the zero point. When , there exists a second non-zero steady-state that reflects well the observed viral behavior [26], and that can be regarded as the analogous of the quasi-stationary state of the exact stochastic SIS model. NIMFA yields an upper bound for the probability of infection of each node, as well as a lower bound for the epidemic threshold [7, 27]. This fact ensures that allows us to determine a safety region for the effective spreading rate, that guarantees the extinction of epidemics in a reasonable time frame. Thus, even though NIMFA is approximated, a design for our optimization problem, based on NIMFA, is always “safe” or “secure”.

In this section, we consider a heterogeneous setting. We include the possibility that the infection rate is link specific, denoting by the infection rate of node towards node . Moreover, each node recovers at rate , so that the curing rate is node specific. Basically, we allow for the epidemics to spread over a graph.

A direct weighted graph (or weighted digraph) is a triple , where the elements of , named arcs (or directed edges), are ordered couples of distinct vertices of , and is a given function; is called the weight of . The matrix , with elements , is the weighted adjacency matrix of . In our framework, and means that node can infect node with rate . Again, self-loops and multiple edges (multiple arcs with the same direction) are not permitted. Hereafter, we shall assume that the directed graph is strongly connected, i.e., for all pairs of nodes , there is a path form to and from to .

As in the homogeneous case, the SIS model with heterogeneous infection and recovery rates is a Markovian process. The time for infected node to infect any susceptible neighbor is an exponential random variable with mean . Also, the time for node to recover is an exponential random variable with mean . A NIMFA model for the heterogeneous setting has been presented first in [28], where a node can infect all neighbors with the same infection rate . Here we include the possibility that the infection rates depend on the connection between two nodes, thus covering a much more general case. The NIMFA governing equation for node in the heterogeneous setting writes as

 (2.1)

Let the vector and let be the matrix defined by when , and ; moreover let be a column vector whose -th component is . Then we can rewrite (2.1) in the following form:

 dP(t)dt=¯¯¯¯AP(t)+F(P). (2.2)

Let

 r(¯¯¯¯A)=max1⩽j⩽NRe(λj(¯¯¯¯A)),

be the stability modulus [29] of , where denotes the real part of the eigenvalues of , . Now, we adopt a result from [29] that lead us to find the epidemic threshold, and to extend the global stability analysis of the homogeneous NIMFA system (see e.g. [20]) to the entirely heterogeneous setting, where each node can potentially infect each of its neighbors with different infection rates. We underline that to use the result in [29], the matrix needs to be irreducible, this is equivalent to say that its associated digraph must be strongly connected.

###### Theorem 2.1

If then is a globally asymptotically stable equilibrium point in for the system (2.2), On the other hand if then there exists a constant solution , such that is globally asymptotically stable in for (2.2).

###### Proof.

See [29, Thm. 3.1].

###### Remark 2.1

Let be an irreducible and non negative matrix, a diagonal matrix with positive entries. Let be the spectrum of the matrix , then the eigenvalue , such that , is real (this follows also from (7.6)).

The result in Theorem  2.1 is crucial for the cost-optimal curing problem described in the next section. In fact, it identifies the value of the epidemic threshold, separating an absorbing phase, where the epidemics will go extinct, from an endemic phase. Thus, this critical threshold is recognized as a key value for treatment strategies against viral infection.

Now, leveraging on the result in Theorem 2.1, we address the problem of suppressing an epidemic spreading, by a cost-optimal distribution of resources within a networked population. Allocating more resources at each node aims to increase its curing rate, that is reflected, e.g., by speeding up its detection capabilities and treatments. We consider that recovery resources have an associated cost that might be different for each node. Thus, let us define a cost function which measures the expenditure in order to distribute curing resources to all nodes. Let be a real, linear and monotonically increasing function with respect to , whose value represents the effort of modifying the recovery rate of node .

This model fits the case of disease treatment plans: policy makers can distribute different amount of resources (e.g. money for medicines, medical and nursery staff, etc,…) in a network of hospitals, or they can design a different health program for different districts, cities, or nations in the case of a timely mass prophylaxis plan. For instance, in the US, pharmaceutical supply caches and production arrangements have been pre-designated. This is done in order to be used for large-scale ongoing prophylaxis and/or vaccination campaigns in case of sudden intentional or natural outbreaks disease [30].

For now, we take into account an arbitrary weighted network. In Sec. Optimal curing policy for epidemic spreading over a community network with heterogeneous population, instead, we shall provide a cost-optimal curing policy for a network with community structure. Hereafter, we consider , with , that we shall call the cost coefficients, for . Thus, the cost function is the cumulative sum over the nodes’ set

 U(\@fontswitch\mathnormalΔ)=N∑i=1ciδi,

where is the curing rate vector.

Now, let us assume that , for all , i.e., the weighted adjacency matrix is symmetric and, consequently, all its eigenvalues are real. Basically, now we are considering undirected graphs with symmetric weights.

Let us define the curing rate matrix, . We remark that, hereafter, we shall indicate with the maximum eigenvalue of . By Theorem 2.1, we know that if , then the epidemics will go extinct. As we have explained in Section Optimal curing policy for epidemic spreading over a community network with heterogeneous population, the critical threshold for the mean-field model is a lower bound of the threshold of the exact Markov model. Thus, the condition corresponds, in the exact stochastic model, to a region where the infectious process dies out exponentially fast for sufficiently large times [24]. We recall that, instead, above the exact threshold the overall-healthy state is reached only after an unrealistically long time. Hence, in order to find a cost-optimal distribution of resources that guarantees the extinction, we seek for the solution of the following problem.

###### Problem 3.1 (Eigenvalue Constraint Formulation)

Find which solves

 minimizeU(\@fontswitch\mathnormalΔ) {subject to:}λ1(A−diag(\@fontswitch\mathnormalΔ))⩽0,\@fontswitch\mathnormalΔ⩾0.

Problem 3.1 can be reformulated as a semidefinite program, that is a convex optimization problem [31]. In fact , where is the -th component of and is the -th element of the standard basis so that . Hereafter, as in [32], the inequality sign in , when is a matrix, means that is positive semidefinite. Thus, we can express the optimization problem with eigenvalue constraint as a semidefinite programming (SDP) problem.

###### Problem 3.2 (Semidefinite Programming Formulation)

Find which solves

 minimizeU(\@fontswitch\mathnormalΔ) {subject to:}diag(\@fontswitch\mathnormalΔ)−A⩾0 {subject to:}\@fontswitch\mathnormalΔ⩾0

The feasibility of the problem is always guaranteed, as showed in the following

###### Theorem 3.3 (Feasibility)

Problem 3.2 is feasible.

###### Proof.

We define and choose , where is the all-one vector of length , consequently, , with identity matrix of order . Then, for any vector , where , for and is an eigenvector basis of , it holds

 wT(A−D)w=wT(N∑i=1λi(A)zivi−lmaxw)⩽(λ1(A)−lmax)||w||2⩽0,

where the last inequality follows since . Hence the chosen vector satisfies the constraint and we can assert that the feasible region is not empty.

Since the problem is feasible there is always an optimal point on the boundary [31] and, by the fundamental result of convex optimization, any locally optimal point of a convex problem is globally optimal [32, Sec. 4.4.2].

Existing results. As introduced in Sec. Optimal curing policy for epidemic spreading over a community network with heterogeneous population, an SDP approach is adopted also in [8] to detect a cost-optimal distribution of protective resources in an arbitrary undirected network. Unlike our approach, they consider that each node can infect all its neighbors with the same infection rate ; moreover they describe the minimization of a decreasing vaccination cost function, which depends on the infection rates, that are allowed to be in a feasible interval. In the second part of the work they propose a greedy approach for the case of all-or-nothing vaccination, i.e., they restrict the infection rate to be in a discrete set, possibly different for each node, , where the two values, are fixed a priori.

With respect to their approach, in our model, each node can potentially infect each of its neighbors with different infection rates, thus we treat a wider scenario. In addition, from the next section onwards, we shall focus, mainly, on a population divided into communities, obtaining a dimensionality reduction of the optimization problem (3.2). Moreover, in Sec. Optimal curing policy for epidemic spreading over a community network with heterogeneous population we propose a bisection algorithm for a two-level optimal curing problem, i.e we consider a two-dimension curing policy, providing that the population is divided in two categories, each of which will benefit from one of the two policies. The two available values of the curing rate are not fixed a priori.

At last, in [9], Preciado et al., leverage on Geometric Programming (GP) techniques for the resource allocation problem, applied to arbitrary weighted directed graphs, hence they do not require the symmetry of the adjacency matrix. However, the drawback of such formulation is that it does not fit for a linear cost function of the type we are considering, which is, anyway, a standard cost function of practical relevance. Thus, in the next section, we show how our formulation of the problem, involving a linear cost function, can be extended to a certain class of not symmetric matrices.

The formulation of the optimization problem (3.2) holds for symmetric weighted adjacency matrix, however we shall show how it can be extended to a certain class of not symmetric matrices that are diagonally symmetrizable. In this case, for a not symmetric matrix , there exists a diagonal matrix such that is symmetric, for similarity their eigenvalues are the same and the semidefinite program formulation can be applied to 111As suggested in  [2], see [33] for a method to check if a matrix is symmetrizable, and, in case, the way to chose the diagonal matrix to achieve symmetry.. Thus, we can include also the case of not symmetric weighted adjacency matrix. A notable example is that of an undirected network where each node can infect all its neighbors with the same infection rate : the weighted adjacency matrix with symmetric and is not symmetric, however it is diagonally symmetrizable, indeed choosing we have , which is symmetric (see, e.g., [8]). Hence, for our problem 3.2, we can request that the matrix is semidefinite positive.

Otherwise, if we have an arbitrary, strongly connected, directed weighted graph and a not symmetrizable -dimensional adjacency matrix , we can apply our formulation to its Hermitian part, , obtaining a suboptimal solution. More precisely, let be the vector of the real part of the eigenvalues of and , both arranged in decreasing order; then it holds [34, Thm. 10.28], meaning that majorizes [34, Sec. 10.1]. Basically, this result suggests that the stability modulus (see Remark 2.1) is smaller than , hence the feasible region of our optimization problem, i.e., where , is a subset of the feasible region for the matrix . Hence, if we solve the problem (3.2), choosing , the value of the cost function obtained is an upper bound of the cost function that would be sufficient to bring at the critical value zero. Thus, we obtain a suboptimal solution, i.e., we will be able to lead the epidemic towards the extinction, but with more effort than it would be sufficient.
Hence, let us consider a diagonally symmetrizable weighted adjacency matrix ; we want to compare the optimal cost function – corresponding to the optimal solution of the problem (3.2)) – with the suboptimal cost function, obtained considering the hermitian part of . Besides, we compute also the cost in the case of a uniform curing rate vector for which the maximum eigenvalue of attains zero. We use a standard solver for semidefinite programs (see [35]). In Fig. 1 a) we consider the cost functions obtained averaging over 300 instances of Erdős-Rényi random graphs with for increasing values of . We take a matrix , where the infection rates are generated as uniform random variables in the interval , and the ’s constants in the interval . We observe that the suboptimal cost function is close to the optimal cost function, the closer the lower the values of . In Fig. 1 b), instead, we fix the value of and we plot the costs as functions of the number of nodes . We can see a growth in the difference between the suboptimal and the optimal cost functions as the number of nodes increase. Ultimately, we obtain always an advantage in the use of suboptimal cost function with respect to the uniform case. In the rest of the paper, we shall consider the case of a network with community structure. We shall show that – in order to find the epidemic threshold for the system (2.2) – it can be employed a matrix with lower dimension than the starting -dimensional adjacency matrix. In turn, this provides a corresponding reduction in the dimension of our optimization problem (see Sec. Optimal curing policy for epidemic spreading over a community network with heterogeneous population).

Hereafter, we shall focus on the case when a contact network is shaped by an existing community network. This framework captures some of the most salient structural inhomogeneities in contact patterns in many applied contexts [36]. There exists an extensive literature on the effect of network community structures on epidemics. A specific community structure may arise due to, for example, geographic separation. Models utilizing this structure are commonly known as “metapopulation” models, where the population is compound of multiple interacting groups, which internally have homogeneous mixing [2] (see, e.g.,  [37, 38]). Such models assume that each community shares a common environment or is defined by a specific relationship. Some of the most common works on metapopulation regard a population divided into households with two level of mixing ([39, 40, 41]). This model typically assumes that contacts, and consequently infections, between nodes in the same group occur at a higher rate than those between nodes in different groups [36]. Thus, groups can be defined, e.g., in terms of spatial proximity, considering that between-group contact rates (and consequently the infection rates) depend in some way on spatial distance, so that, each individual can be theoretically infected by each of the other individuals. However, models where infection can only be transmitted by nodes directly connected by an edge, may provide a more realistic approach to the study of the evolution of the epidemics. In turn, an important challenge is how to consider a realistic underlying structure and appropriately incorporate the influence of the network topology on the dynamics of epidemic [36, 42, 43, 44, 45].

In [20, 46] the authors analyze the dynamics of an epidemics on networks that are partitioned into local communities, through the first-order mean-field approximation discussed in Section Optimal curing policy for epidemic spreading over a community network with heterogeneous population. The investigation was based on the graph-theoretical notion of            equi [14, 15, 16]. Specifically, for an undirected graph, let be a partition of the node set , i.e., a sequence of mutually disjoint nonempty subsets of , called cells, whose union is , that we assume given a priori; is called equitable if the subgraph of induced by is regular for all ’s. Furthermore, for any two subgraphs and , whenever there exists at least one connection between nodes in the first and second subgraph, then each node in is connected with the same number of nodes in . In [20, 46] two-level infection rates have been considered: an intra-community infection rate and a lower inter-community infection rate. In the network structures hereafter, we generalize the model to more than two levels of infectiousness. We observe that usually a community is defined as a set of network nodes joined together in tightly-knit groups whereas among such group connections are looser [47]. Leveraging on the definition of equitable partition, instead, we can also consider that connections between nodes belonging to the same community can be, eventually, less dense than connections with other communities. Thus, the definition of community acquires a broader sense.

Networks with an equitable partition of the node set can describe models consisting of multiple smaller sub-populations such as, e.g., households, workplaces, or classes in a school, when the internal structure of each community is represented by a complete graph (members of a small community usually know each other) and all the nodes of adjacent communities are mutually linked (all member of adjacent communities may potentially come into contact). Equitable partitions can be observed also in the architecture of some computer networks where clusters of clients connect to single routers, whereas the router network has a connectivity structure with the nodal degree constrained by the number of ports (see as examples Fig. 2b). Equitable partitions appear also in the study of synchrony and pattern formation in coupled cell networks [48, 49] where they are named âbalancedâ partitions. Equitable partitions have been used also to analyze the controllability of multi-agent systems, for the case of a multi-leader setting [50], and for the leader-selection controllability problem, in characterizing the set of nodes from which a given networked control system (NCS) is controllable/uncontrollable [51]. These works show interesting realistic scenarios for the use of equitable partitions.

In particular, since the size of some real networks might pose limitations in our ability to investigate their spectral properties, we can leverage on the structural regularity of network with equitable partition, to reduce the dimensionality of our optimization problem (3.2).

Next, we define equitable partitions for the case of a directed weighted networks, extending the analysis in [20] to this framework. With a little abuse of notation, hereafter we shall refer to a partition of a network, to indicate the partition of its node set.

The definition of equitable partitions can be extended to weighted directed graphs, based on [16, Def. 8.24]. That definition applies to oriented weighted graphs [16, Def. A.1]: we prefer to allow for a pair of symmetric oriented edges in order to cover naturally unoriented graphs.

###### Definition 4.1

Let be a weighted directed graph. The partition of the node set is called or if for all , there are

 cinij∈Rs.t.∑w∈Vjρ((v,w))=cinij,for allv∈Vi,

or

 coutij∈Rs.t.∑w∈Vjρ((w,v))=coutij,for allv∈Vi,

respectively. The partition is called if it is both inward and outward equitable, hence for all , there are

 cij∈Rs.t.∑w∈Vjρ((v,w))+ρ((w,v))=cij,for allv∈Vi.

We shall identify the set of all nodes in with the -th of the whole population.

###### Remark 4.1

Let be the number of elements of , . If the partition of the node set of a weighted di-graph is equitable, then for all ,

 kicoutij=kjcinji, (4.1)

An equitable partition generates the , which is a , directed and weighted, with cells as vertices. For the sake of explanation, in the following, we will identify with the simple graph having the same vertex set, i.e. composed by the cells, where an edge exists between two cells, if at least one exists in the original multigraph.

For the purpose of modeling, nodes of the quotient graph can represent communities, e.g., villages, cities or countries. Link weights in the quotient graph, in turn, provide the strength of the contacts between such communities. In particular, the weight of a link may be (a non-negative) function of the number of people traveling per day between two countries; in fact, the frequency of contacts between them correlates with the propensity of a disease to spread between nodes.

Related to the quotient graph, there exists a quotient matrix, that contains the relevant structural information of the networks. Thus, let us consider the matrix , where

 siv={1√|Vi| v∈Vi0otherwise,

from which it follows that . Now let us consider the transpose of the adjacency matrix of the weighted directed graph , that is

 AT=¯¯¯¯A+D, (4.2)

where, we remember, is the matrix in (2.2) and is the curing rate matrix. Then, the transpose of the of (with respect to the given partition) is

 QT:=SATST.

We can write the following explicit expression for :

 QT=diag(coutii)+(kicoutij√kikj)i,j=1,…,n (4.3)

By (4.1), we can write the matrix as

 QT=diag(coutii)+(√coutijcinji)i,j=1,…,n (4.4)
###### Remark 4.2

We observe that matrix in (4.4) might not be symmetric, whereas in the case of undirected graphs it is always symmetric (see, e.g., [15, 20]). Even though we have represented the most general definition of an equitable partition simpler situations can be represented. E.g., nodes of the same community may infect all nodes in another community with the same rate.

When considering a population partitioned into communities, it may be appropriate to take into account the case where all nodes of a tagged community have the same recovery rate , . In turn, such rate may differ from one community to the other. We remember that is the total number of nodes in the network, whereas is the number of communities.

###### Definition 4.2

Let us introduce the vector of nonzero curing rates , that we shall call the reduced curing rate vector and , the reduced curing rate matrix.

Thus, we have the curing rates vector , with components for all and .

In appendix Optimal curing policy for epidemic spreading over a community network with heterogeneous population we shall discuss when and how it is possible to reduce the original system (2.2) to a system of differential equations, through the matrix . Since for our optimization problem the parameter of interest is the epidemic threshold, in this section we limit our self to results related to this critical value.

###### Lemma 4.1

Let be an equitable partition. Let and be weighted matrices as in (4.2) and (4.4), respectively. Then, it holds:

i) .
ii) For all and all

 (QT−¯¯¯¯¯D)x=λxif and only if(AT−D)STx=λSTx.

Now let us consider the system of differential equations (2.2). It is possible to extend [20, Thm. 4.1] to the case of directed graphs. Following that result, if we assume that at time the infection probability is equal for all nodes in the same community (while it may differ from one community to the other), the number of equations in (2.2) can be reduced by using the transpose of the quotient matrix . Hence, the reduced dynamical system writes

 d¯¯¯pj(t)dt =(1−¯¯¯pj(t))n∑m=1m≠jcoutjm¯¯¯pm(t)+coutjj(1−¯¯¯pj(t))¯¯¯pj(t)−¯¯¯¯δj¯¯¯pj(t),j=1,…,n (4.5)

where is the infection probability of a node in the community . We can prove that, in the case of a graph whose node set has an equitable partition, and , the critical threshold for (2.2), applying Thm. 2.1, can be determined directly considering the reduced system (7.1).

###### Proposition 4.3

The elements of the curing rates vector , that determine the critical threshold of (2.2), are identified by the elements of , in such a way that for all , , for which

 r(QT−¯¯¯¯¯D)=0, (4.6)

where is the stability modulus.

Since the quotient matrix and the adjacency matrix have the same stability modulus (and so their transposed do), a computational advantage can be obtained in the calculation of the critical threshold of the system (2.2). This result is very relevant for our optimization problem. Indeed, in the case of a network with equitable partition, we can use a lower dimensional matrix to compute the epidemic threshold.

In this section, we consider a heterogeneous curing control per community. First, we assume that all nodes in community infect all nodes in community with the same infection rate, , and that , . In this case, the graph is undirected and the weights are symmetric, thus the quotient matrix is symmetric and has real eigenvalues. Now let us consider the reduced curing rate vector , the cost function writes

 U(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\@fontswitch\mathnormalΔ)=n∑j=1cjkj¯¯¯δj.

Thus, is the cost for curing all elements of each community at rate , where , . We seek for the solution of the following

###### Problem 5.1 (Eigenvalue Constraint Formulation)

Find which solves

 minimizeU(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\@fontswitch\mathnormalΔ) {subject to:}λ1(Q−diag(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\@fontswitch\mathnormalΔ))⩽0,¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\@fontswitch\mathnormalΔ⩾0

which also writes

###### Problem 5.2 (Semidefinite Programming Formulation)

Find which solves

 minimizeU(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\@fontswitch\mathnormalΔ) {subject to:}diag(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\@fontswitch\mathnormalΔ)−Q⩾0 {subject to:}¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯\@fontswitch\mathnormalΔ⩾0

Thm. 3.3 guarantees the feasibility of the problem. The general case of equitable partitions introduced in Sec. Optimal curing policy for epidemic spreading over a community network with heterogeneous population, may not lead to a symmetric quotient matrix . However, we may consider suboptimal solutions – as explained in Sec. Optimal curing policy for epidemic spreading over a community network with heterogeneous population – obtained by applying the formulation of our optimization problem to the hermitian part of .

In the next section, we consider a simpler version of Problem  5.2 and we design a more efficient algorithm with respect to the SDP program.

The state of the art for SDP solvers such as, e.g., the SDPT3 solver used for our numerical computation, provide solutions for moderate size graphs. Actually, the best known bound for the complexity of an -solution attained with an interior point method is , where represents the accuracy [52]. The problem can be solved more efficiently when we face a two-level optimal curing problem, for which we shall provide an algorithm that yields an -approximation of the optimal solution, with a complexity equal to (see Thm. 5.7). Precisely, we consider only two possible levels of the nodes local curing rates, let us say and , that are not fixed a priori. This situation fits well, e.g., in the case of a network where communities are of “two types”. Communities of the first type are eligible for curing rate , whereas communities of the second type are eligible for curing rate . For convenience, we define the former, communities, and the latter, communities.

Such kind of configuration is suitable for a network that is, e.g., bipartite (where each node, e.g., represents a full-meshed community), or for an interconnected stars network, i.e., a network obtained by interconnecting star graphs by linking stars’ central nodes (see Fig 2b). Let us note that the Barabási-Albert graph model [53], that captures the power-law degree distribution often seen (or approximately seen) in real-world networks, can be regarded as a set of hubs with star graph features [54]. Bipartite networks, instead, can be used to understand the spreading of sexually transmitted diseases, in which the population is naturally divided into males and females and the disease can only be transmitted between nodes of different kinds. Bipartite networks can also represent the spreading of diseases in hospitals, in which one type of node accounts for (isolated) patients and the other type for caregivers, or some vector-borne diseases, such as malaria, in which the transmission can only take place between the vectors and the hosts [1].

Thus, let us consider the following partition of the node set, and . We assume that the node set partition is equitable. Let us introduce the curing matrix and define

 I0m=[Im000],I1m′=[000Im′]

where is the identity matrix of order . Then, we can write the semidefinite programming for the two-level curing rates, shortly the curing problem, as follows:

###### Problem 5.3 (Semidefinite Programming 2D Formulation)

Find which solves

 minimizeU(\@fontswitch\mathnormalΔ2) {subject to:}δ0I0m+δ1I1m′−Q⩾0 {subject to:}\@fontswitch\mathnormalΔ2⩾0

The cost function is

 U(\@fontswitch\mathnormalΔ2)=∑Vj∈π0kjf0(δ0)+∑Vz∈π1kzf1(δ1),

where , and , with , represent the effort to modify the recovery rate for nodes belonging to , and , respectively.

In Section Optimal curing policy for epidemic spreading over a community network with heterogeneous population of the Appendix, we shall provide some simple examples for the optimal solution of the Problem (5.3).

In the design of our algorithmic solution, we have leveraged on some basic properties of the curing problem. In order to do so, we need a few basic facts recalled next.

###### Proposition 5.4

Let be an symmetric, irreducible and non negative matrix and let :
i. if for some , then ;
ii. The function is continuous;
iii. is strictly decreasing in , .

Let us denote by the feasibility region of Prob. 5.3: same argument of Thm. 3.3 let us state that it is non empty; furthermore, the problem structure guarantees to be convex [32]. We indicate and the standard projections of onto the -axis and the -axis, respectively.

###### Lemma 5.1 (Monotonicity)

Let be the function that associates to the value such that . Then, is decreasing.

###### Proof.

First, let us show that is a well defined function over . Let , because of feasibility, there exists , where , such that . Furthermore, it holds by of Prop. 5.4. Because of and in Prop. 5.4 we know that is a continuous strictly decreasing function of over , so there exists one and only one value satisfying the definition of .

Let and assume that , for some , i.e, that is not decreasing. From the definition of there exists . Hence, we can write

 wT(Q−diag(δ01m,ϕ(δ0)1m′))w=wTdiag(z1m,ζ1m′)w+wT(Q−diag((δ0+z)1m,ϕ(δ0+z)1m′))w =wTdiag((z1m,ζ1m′))w>0,

where the strict inequality holds because . Since , this means that must be semidefinite negative and we have a contradiction.

We prove next that the search for the optimal solution can be restricted to a compact subset of .

###### Theorem 5.5 (Compact search set)

There exist two pairs and such that a solution of Prob. 5.3 belongs to a compact subset .

###### Proof.

Let us define and , then we write , with as in Thm. 3.3, and . Let us denote , by Thm. 3.3, , hence and, by defining set , it follows that ; is closed as intersection of closed sets.

Now, feasibility conditions of Prob. 5.3 require matrix to be semidefinite negative. We define : we have since and by of Prop. 5.4. By assertion in Prop. 5.4, is a continuous function. Hence, there exists such that , and since is decreasing . We can repeat the same reasoning by inverting the role of and defining . Hence, we can assert that exists such that and .

Finally, by letting , the points and belong to , i.e., they belong to , so , and consequently, being closed, it is also compact.

###### Remark 5.1

Thm. 5.5 allows us to identify an interval of the values of and where we can restrict the search of . Since