TY - JOUR
AU - Van Mieghem, Piet,
AB - 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 analyse 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. 1. Introduction 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 behaviours, diffusion of information, computer viruses, etc., indeed, although the basic mechanisms of these phenomena can be different, often their dynamical behaviour 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–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, for example 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 Sahneh and Scoglio [6] take advantage of the $$N$$-intertwined approximation [7] to analyse 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 Section 3, we analyse 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 neighbours. 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. 1.1 Outline and main results 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, for example 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 per community (i.e. per hospital, school, village, or city, etc, $${\ldots}$$), 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 equitable partition [14–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–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 4). The article is organized as follows. In Section 2, first, we review some background concepts for epidemic processes on networks in the homogeneous setting and the related mean-field approximation adopted in the article. 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 Section 3, 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 Section 4, 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, for example when the population can be divided in two categories (younger and elders, male and female, etc, $${\ldots}$$). This kind of situation fits well certain networks with equitable partition, such as, for example bipartite graphs and interconnected star networks. In this case, we provide a scalable bisection algorithm, that yields an $$\varepsilon$$-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 Appendix. 2. The epidemic network model 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 $$G=(V,E)$$, with edge set $$E$$ and node (vertex) set $$V$$. The order of $$G$$, denoted by $$N$$, is the cardinality of $$V$$. The edge set of $$G$$ consists of unordered pairs $$\left\{i,j\right\}$$, with $$i,j \in V$$, and $$i \neq j$$. Connectivity of graph $$G$$ is conveniently expressed by the symmetric $$N \times N$$ adjacency matrix $$A$$. The viral state of a node $$i$$, at time $$t$$, is described by a Bernoulli random variable $$X_i(t)$$, where we set $$X_i(t) = 0$$, if $$i$$ is healthy and $$X_i(t) = 1$$, if $$i$$ is infected. Every node at time $$t$$ is either infected with probability $$p_i(t) = {\mathbb P}(X_i(t) = 1)$$ or healthy (but susceptible) with probability $$1 - p_i(t)={\mathbb P}(X_i(t) = 0) $$. In the homogeneous setting, the recovery process is a Poisson process with rate $$\delta$$, and the infection process is a per link Poisson process where the infection rate between an healthy and an infected node is $$\beta$$. All the infection and recovery processes are independent. The SIS process, developing a graph with $$N$$ nodes, can be modelled as a continuous-time Markov process with $$2^N$$ states, covering all possible combinations in which $$N$$ 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 $$2^N$$ linear differential equations by $$N$$ non-linear differential equations representing the time-change of the infection probability of each node. Epidemic threshold: For a network with finite order $$N$$, 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 behaviour: indeed, there is a critical value $$\tau_c$$ of the effective spreading rate $$\tau= \beta/\delta$$, whereby if $$\tau$$ is lesser than $$\tau_c$$, the initial infection dies out quickly. Conversely, for $$\tau$$ larger than $$\tau_c$$, the infection spreading can last very long in any sufficiently large network [22–24]. The regime of persistent infection ($$\tau > \tau_c$$ ), 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 $$(N \simeq 100)$$ and when $$\tau > \tau_c$$, 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 $$\tau^{(1)}_c =\frac{1}{\lambda_{1}(A)}$$, where $$\lambda_1(A)$$ is the spectral radius of the adjacency matrix $$A$$, (see [7, 25]). When $$\tau \leq \tau_c^{(1)}$$ the only equilibrium of the NIMFA system is the zero point. When $$\tau > \tau^{(1)}_c$$, there exists a second non-zero steady-state that reflects well the observed viral behaviour [25], 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, 26]. This fact ensures that $$\tau_c^{(1)}$$ allows us to determine a safety region $$\left\{\tau \leq \tau^{(1)}_c\right\}$$ 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’. 2.1 Heterogeneous SIS mean-field model In this section, we consider a heterogeneous setting. We include the possibility that the infection rate is link specific, denoting by $$\beta_{ij}$$ the infection rate of node $$j$$ towards node $$i$$. Moreover, each node $$i$$ recovers at rate $$\delta_i$$, so that the curing rate is node specific. Basically, we allow for the epidemics to spread over a directed weighted graph. A direct weighted graph (or weighted digraph) is a triple $$G=(V,E, \rho)$$, where the elements of $$E$$, named arcs (or directed edges), are ordered couples $$e=(i,j)$$ of distinct vertices of $$V$$, and $$\rho:E \rightarrow (0,\infty)$$ is a given function; $$\rho(e)$$ is called the weight of $$e$$. The matrix $$A=(a_{ij})$$, with elements $$a_{ij}=\rho(i,j)=\beta_{ji}$$, is the weighted adjacency matrix of $$G$$. In our framework, $$e=(i,j)\in E$$ and $$\rho(e)=\beta_{ji}$$ means that node $$i$$ can infect node $$j$$ with rate $$\beta_{ji}$$. 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, that is for all pairs of nodes $$i,j \in V$$, there is a path form $$i$$ to $$j$$ and from $$j$$ to $$i$$. As in the homogeneous case, the SIS model with heterogeneous infection and recovery rates is a Markovian process. The time for infected node $$j$$ to infect any susceptible neighbour $$i$$ is an exponential random variable with mean $$\beta_{ij}^{-1}$$. Also, the time for node $$j$$ to recover is an exponential random variable with mean $$\delta_j^{-1}$$. A NIMFA model for the heterogeneous setting has been presented first in [27], where a node $$i$$ can infect all neighbours with the same infection rate $$\beta_i$$. 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 $$i$$ in the heterogeneous setting writes as \begin{equation} \begin{aligned}\label{het1} \hskip-1mm \frac{d p_i(t)}{dt}= \sum_{j=1}^N \beta_{ij} p_j(t) - \sum_{j=1}^N \beta_{ij} p_i(t) p_j(t) -\delta_i p_i(t), \; \hskip15mm i= 1,\ldots,N . \end{aligned} \end{equation} (2.1) Let the vector $$P(t)=(p_1(t), \ldots, p_N(t))^{\top}$$ and let $$\overline{A}=(\overline{a}_{ij})$$ be the matrix defined by $$\overline{a}_{ij}=\beta_{ij}$$ when $$i \neq j$$, and $$\overline{a}_{ii}= - \delta_i$$; moreover, let $$F(P)$$ be a column vector whose $$i$$th component is $$- \sum_{j=1}^N \beta_{ij} p_i(t) p_j(t)$$. Then, we can rewrite (2.1) in the following form: \begin{equation}\label{het} \frac{d P(t)}{dt}= \overline{A} P(t) + F(P). \end{equation} (2.2) Let \[ r(\overline A)= \max_{1 \leq j \leq N} Re(\lambda_j(\overline{A})), \] be the stability modulus [28] of $$\overline A$$, where $$Re (\lambda_j(\overline{A}))$$ denotes the real part of the eigenvalues of $$\overline A$$, $$j=1, \ldots, N$$. Now, we adopt a result from [28] 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 neighbours with different infection rates. We underline that to use the result in [28], the matrix $$\overline{A}$$ needs to be irreducible, this is equivalent to say that its associated digraph must be strongly connected. Theorem 2.1 If $$r(\overline{A}) \leq 0$$ then $$P=0$$ is a globally asymptotically stable equilibrium point in $$I_N =[0,1]^N$$ for the system (2.2), On the other hand if $$r(\overline{A}) > 0$$ then there exists a constant solution $$P^{\infty} \in I_N -\left\{0\right\}$$, such that $$P^{\infty}$$ is globally asymptotically stable in $$I_N -\left\{0\right\}$$ for (2.2). Proof. See [28, Theorem 3.1]. □ Remark 2.1 Let $$A$$ be an $$N \times N$$ irreducible and non-negative matrix, $$D$$ a diagonal matrix with positive entries. Let $$\sigma(A-D)$$ be the spectrum of the matrix $$A-D$$, then the eigenvalue $$\lambda \in \sigma(A-D)$$, such that $$Re(\lambda)=r(A-D)$$, is real (this follows also from (A.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. 3. Optimal curing policies for arbitrary weighted networks 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, for example 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 $$f_i(\delta_i)$$ be a real, linear and monotonically increasing function with respect to $$\delta_i$$, whose value represents the effort of modifying the recovery rate of node $$i$$. 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, $${\ldots}$$) 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 USA, 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 [29]. For now, we take into account an arbitrary weighted network. In Section 4, instead, we shall provide a cost-optimal curing policy for a network with community structure. Hereafter, we consider $$f(\delta_i)=c_i \delta_i$$, with $$c_i>0$$, that we shall call the cost coefficients, for $$i=1,\ldots, N$$. Thus, the cost function is the cumulative sum over the nodes’ set \begin{equation*}\label{eq:cost} U(\Delta)=\sum_{i=1}^N c_i \delta_i, \end{equation*} where $$\Delta=(\delta_1, \ldots, \delta_N)$$ is the curing rate vector. 3.1 Undirected graphs with symmetric weights Now, let us assume that $$\beta_{ij}=\beta_{ji}$$, for all $$i,j=1, \ldots, N$$, that is the weighted adjacency matrix $$A=(\beta_{ij})$$ is symmetric and, consequently, all its eigenvalues are real. Basically, now we are considering undirected graphs with symmetric weights. Let us define the $$N \times N$$ curing rate matrix, $$D=\mathop{\rm diag}(\Delta)$$. We remark that, hereafter, we shall indicate with $$\lambda_1(A)$$ the maximum eigenvalue of $$A$$. By Theorem 2.1, we know that if $$\lambda_1\!\big ( A - \mathop{\rm diag}{( \Delta)} \big ) \leq 0$$, then the epidemics will go extinct. As we have explained in Section 2, the critical threshold for the mean-field model is a lower bound of the threshold of the exact Markov model. Thus, the condition $$\lambda_1\!\big ( A - \mathop{\rm diag}{( \Delta)} \big ) \leq 0$$ 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 $$\ \Delta \geq 0$$ which solves \begin{eqnarray*} && {\rm{minimize}} \qquad\qquad U(\Delta)\\ &&\mbox{subject to:} \quad\;\; \lambda_1\!\big ( A - \mathop{\rm diag}{( \Delta)} \big ) \leq 0, \quad\;\; \Delta \geq 0 \nonumber. \end{eqnarray*} Problem 3.1 can be reformulated as a semidefinite program, that is a convex optimization problem [30]. In fact $$\mathop{\rm diag}(\Delta)=\sum_{i=1}^{N} \Delta_i \mathop{\rm diag}(\mathbf{e}_i)$$, where $$\Delta_i$$ is the $$i$$th component of $$\Delta$$ and $$\mathbf{e}_i$$ is the $$i$$th element of the standard basis so that $$\mathop{\rm diag}(\mathbf{e}_i)\geq 0$$. Hereafter, as in [31], the inequality sign in $$M \geq 0$$, when $$M$$ is a matrix, means that $$M$$ is positive semidefinite. Thus, we can express the optimization problem with eigenvalue constraint as a SDP problem. Problem 3.2 (SDP formulation) Find $$\Delta$$ which solves \[\begin{align} & \text{minimize}\qquad \qquad \quad U(\Delta ) \\ & \text{subject to:}\quad \quad \text{diag}\,\,(\Delta )\,\,\,-A\,\,\,\,\ge \,\,\,0 \\ & \quad \ \ \quad \quad \quad \quad \quad \Delta \,\,\,\ge \,\,0. \\ \end{align}\] 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 $$l_{\max}:=\max_i\sum_j a_{ij}$$ and choose $$\Delta= l_{\max} \mathbf{1}_{N}$$, where $$\mathbf{1}_{N}$$ is the all-one vector of length $$N$$, consequently, $$D=l_{\max} I_{N}$$, with $$I_{N}$$ identity matrix of order $$N$$. Then, for any vector $$w=\sum_{i=1}^N z_i v_i$$, where $$z_i \in \mathbb R$$, for $$i=1 \ldots N$$ and $$\left\{v_1, \ldots, v_N\right\}$$ is an eigenvector basis of $$A$$, it holds \begin{eqnarray*} && w^{\top}(A- D)w = w^{\top}\left( \sum_{i=1}^N \lambda_i(A)z_iv_i - l_{\max} w\right) \leq ( \lambda_1(A) - l_{\max})||w||^2 \leq 0, \end{eqnarray*} where the last inequality follows since $$\lambda_1(A) \leq \max_i\sum_j a_{ij}$$. 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 [30] and, by the fundamental result of convex optimization, any locally optimal point of a convex problem is globally optimal [31, Section 4.4.2]. Existing results: As introduced in Section 1, 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 $$i$$ can infect all its neighbours with the same infection rate $$\beta_i$$; 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, that is they restrict the infection rate to be in a discrete set, possibly different for each node, $$\beta_i \in \left\{\underline{\beta_i}, \overline{\beta_i}\right\}$$, where the two values, are fixed a priori. With respect to their approach, in our model, each node can potentially infect each of its neighbours 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 Section 5.3, we propose a bisection algorithm for a two-level optimal curing problem, i.e. we consider a two-dimensional 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. 3.2 Extension to directed weighted networks 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 $$A$$, there exists a diagonal matrix $$G$$ such that $$G^{-1}A G$$ is symmetric, for similarity their eigenvalues are the same and the semidefinite program formulation can be applied to $$G^{-1}A G$$1. 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 $$i$$ can infect all its neighbours with the same infection rate $$\beta_i$$: the weighted adjacency matrix $$BA$$ with $$A$$ symmetric and $$B=\mathop{\rm diag}(\beta_i)$$ is not symmetric; however, it is diagonally symmetrizable, indeed choosing $$G=B^{1/2}$$ we have $$B^{-1/2}(BA) B^{1/2}= B^{1/2} A B^{1/2}$$, which is symmetric (see, e.g. [8]). Hence, for our Problem 3.2, we can request that the matrix $$B^{1/2} A B^{1/2}$$ is semidefinite positive. Otherwise, if we have an arbitrary, strongly connected, directed weighted graph and a not symmetrizable $$N$$-dimensional adjacency matrix $$A$$, we can apply our formulation to its Hermitian part, $${\mathcal H}=(A+A^{\top})/2$$, obtaining a suboptimal solution. More precisely, let $$Re\lambda(A)=\left(Re \lambda_1(A), \ldots, Re \lambda_N(A) \right)$$ be the vector of the real part of the eigenvalues of $$A$$ and $$\lambda({\mathcal H} (A))=(\lambda_1({\mathcal H} (A)), \ldots, \lambda_N({\mathcal H} (A)))$$, both arranged in decreasing order; then it holds $$Re\lambda(A) \prec \lambda({\mathcal H} (A))$$ [33, Theorem 10.28], meaning that $$\lambda({\mathcal H} (A))$$majorizes$$Re\lambda(A)$$ [33, Section 10.1]. Basically, this result suggests that the stability modulus $$\lambda_1(A-D)$$ (see Remark 2.1) is smaller than $$\lambda_1({\mathcal H}(A)-D)$$, hence the feasible region of our optimization problem, that is where $$\lambda_1({\mathcal H}(A)-D) \leq 0 $$, is a subset of the feasible region for the matrix $$A$$. Hence, if we solve the problem (3.2), choosing $${\mathcal H}(A)$$, the value of the cost function obtained is an upper bound of the cost function that would be sufficient to bring $$\lambda_1(A-D)$$ at the critical value zero. Thus, we obtain a suboptimal solution, that is 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 $$BA$$; 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 $$BA$$. Besides, we compute also the cost in the case of a uniform curing rate vector for which the maximum eigenvalue of $$B^{1/2}AB^{1/2}$$ attains zero. We use a standard solver for semidefinite programs (see [34]). In Fig. 1(a), we consider the cost functions obtained averaging over 300 instances of Erdős–Rényi random graphs with $$N=100$$ for increasing values of $$p$$. We take a matrix $$B=\mathop{\rm diag}(\beta_i)$$, where the infection rates are generated as uniform random variables in the interval $$(0.1,3)$$, and the $$c_i$$’s constants in the interval $$(0.5,5)$$. We observe that the suboptimal cost function is close to the optimal cost function, the closer the lower the values of $$p$$. In Fig. 1(b), instead, we fix the value of $$p=0.3$$ and we plot the costs as functions of the number of nodes $$N$$. 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. Fig. 1. View largeDownload slide Extension to directed weighted networks: comparison between optimal, suboptimal and uniform solutions. (a) the cost values have been obtained averaging over 300 instances of Erdős–Rényi sample graphs with $$N=100$$, for different $$p=0.1,0.3, 0.5,0.7, 0.9$$ and (b) average over $$300$$ Erdős–Rényi sample graphs with $$N=100,200,400,600,800,1000$$ for $$p=0.3$$; $$0.95$$ confidence intervals are superimposed. Fig. 1. View largeDownload slide Extension to directed weighted networks: comparison between optimal, suboptimal and uniform solutions. (a) the cost values have been obtained averaging over 300 instances of Erdős–Rényi sample graphs with $$N=100$$, for different $$p=0.1,0.3, 0.5,0.7, 0.9$$ and (b) average over $$300$$ Erdős–Rényi sample graphs with $$N=100,200,400,600,800,1000$$ for $$p=0.3$$; $$0.95$$ confidence intervals are superimposed. In the rest of the article, 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 $$N$$-dimensional adjacency matrix. In turn, this provides a corresponding reduction in the dimension of our optimization problem (see Section 5). 4. Community networks 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 [35]. 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. [36, 37]). 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 [38–40]. 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 [35]. Thus, groups can be defined, for example 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 [35, 41–44]. In [20, 45] the authors analyse the dynamics of an epidemics on networks that are partitioned into local communities, through the first-order mean-field approximation discussed in Section 2. The investigation was based on the graph-theoretical notion of equitable partition [14–16]. Specifically, for an undirected graph, let $$\pi=\left\{V_1,\,{\ldots}\,,V_n\right\}$$ be a partition of the node set $$V$$, that is a sequence of mutually disjoint non-empty subsets of $$V$$, called cells, whose union is $$V$$, that we assume given a priori; $$\pi$$ is called equitable if the subgraph $$G_i$$ of $$G(V,E)$$ induced by $$V_i$$ is regular for all $$i$$’s. Furthermore, for any two subgraphs $$G_i$$ and $$G_j$$, whenever there exists at least one connection between nodes in the first and second subgraph, then each node in $$G_i$$ is connected with the same number of nodes in $$G_j$$. In [20, 45], 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 [46]. 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, for example 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 [47, 48] where they are named balanced partitions. Equitable partitions have been used also to analyse the controllability of multi-agent systems, for the case of a multi-leader setting [49], and for the leader-selection controllability problem, in characterizing the set of nodes from which a given networked control system (NCS) is controllable/uncontrollable [50]. These works show interesting realistic scenarios for the use of equitable partitions. Fig. 2. View largeDownload slide A sample graphs with equitable partition. (a) $$V=\{V_1, V_2, V_3, V_4\}$$ and (b) interconnected star networks: $$V=\{V_1^0, V_2^0, V_3^0, V_4^0, V_5^0, V_1, V_2, V_3, V_4, V_5\}.$$ Fig. 2. View largeDownload slide A sample graphs with equitable partition. (a) $$V=\{V_1, V_2, V_3, V_4\}$$ and (b) interconnected star networks: $$V=\{V_1^0, V_2^0, V_3^0, V_4^0, V_5^0, V_1, V_2, V_3, V_4, V_5\}.$$ 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. 4.1 Equitable partitions for weighted directed networks The definition of equitable partitions can be extended to weighted directed graphs, based on [16, Definition 8.24]. That definition applies to oriented weighted graphs [16, Definition A.1]: we prefer to allow for a pair of symmetric oriented edges in order to cover naturally unoriented graphs. Definition 4.1 Let $$G=(V,E,\rho)$$ be a weighted directed graph. The partition $$\pi=\left\{V_1,\,{\ldots}\,,V_n\right\}$$ of the node set $$V$$ is called inward equitable or outward equitable if for all $$i,j \in \left\{1, \dots ,n \right\}$$, there are \begin{align*} c^{in}_{ij} \in \mathbb{R} \quad \text{s.t.} \quad \sum_{w \in V_j} \rho((v,w)) = c^{in}_{ij}, \quad \hskip5mm \text{for all} \quad v \in V_i, \end{align*} or \begin{align*} c^{\rm out}_{ij} \in \mathbb{R} \quad \text{s.t.} \quad \sum_{w \in V_j} \rho((w,v)) = c^{\rm out}_{ij}, \hskip8mm \text{for all} \quad v \in V_i, \end{align*} respectively. The partition is called equitable if it is both inward and outward equitable, hence for all $$i,j \in \left\{1, \dots ,n \right\}$$, there are \begin{align*} c_{ij} \in \mathbb{R} \quad \text{s.t.} \quad \sum_{w \in V_j} \rho((v,w)) + \rho((w,v))= c_{ij}, \hskip8mm \text{for all} \quad v \in V_i. \end{align*} We shall identify the set of all nodes in $$V_i$$ with the $$i$$-th community of the whole population. Remark 4.1 Let $$k_i$$ be the number of elements of $$V_i$$, $$i=1, \ldots, n$$. If the partition of the node set of a weighted di-graph is equitable, then for all $$i,j \in \left\{1, \dots ,n \right\}$$, \begin{equation}\label{out} k_i c^{\rm out}_{ij}=k_j c^{in}_{ji}, \end{equation} (4.1) An equitable partition generates the quotient graph$$G/\pi$$, which is a multigraph, directed and weighted, with cells as vertices. For the sake of explanation, in the following, we will identify $$G/\pi$$ with the simple graph having the same vertex set, that is composed by the cells, where an edge exists between two cells, if at least one exists in the original multigraph. For the purpose of modelling, nodes of the quotient graph can represent communities, for example 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 travelling 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 $$n \times N$$ matrix $$S=(s_{iv})$$, where \begin{equation*} s_{iv}=\begin{cases} \frac{1}{\sqrt{|V_i|}} & v\,\in V_i\\ 0 & \text{otherwise}, \end{cases} \end{equation*} from which it follows that $$SS^{\top}=I_n$$. Now let us consider the transpose of the adjacency matrix of the weighted directed graph $$G$$, that is \begin{equation}\label{matb4} A^{\top}=\overline{A}+D, \end{equation} (4.2) where, we remember, $$\overline{A}$$ is the matrix in (2.2) and $$D=\mathop{\rm diag}(\Delta)$$ is the curing rate matrix. Then, the transpose of the quotient matrix of $$G$$ (with respect to the given partition) is \[ Q^{\top}:=SA^{\top}S^{\top}. \] We can write the following explicit expression for $$Q^{\top}$$: \begin{equation} Q^{\top}=\mathop{\rm diag}(c^{\rm out}_{ii})+ \left(\frac{k_i c^{\rm out}_{ij}}{\sqrt{k_ik_j}}\right)_{i,j=1, \ldots, n} \end{equation} (4.3) By (4.1), we can write the matrix $$Q^{\top}$$ as \begin{equation}\label{eq:q4} Q^{\top}=\mathop{\rm diag}(c^{\rm out}_{ii})+ \left(\sqrt{c^{\rm out}_{ij} c^{in}_{ji}}\right)_{i,j=1, \ldots, n} \end{equation} (4.4) Remark 4.2 We observe that matrix $$Q$$ 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. For example, 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 $$j$$ have the same recovery rate $$\delta_j$$, $$j=1, \ldots ,n$$. In turn, such rate may differ from one community to the other. We remember that $$N$$ is the total number of nodes in the network, whereas $$n$$ is the number of communities. Definition 4.2 Let us introduce the $$1 \times n$$ vector of non-zero curing rates $$\overline{\Delta}=(\overline{\delta}_1, \,{\ldots}\, , \overline{\delta}_n)$$, that we shall call the reduced curing rate vector and $$\overline{D}=\mathop{\rm diag}(\overline{\Delta})$$, the reduced curing rate matrix. Thus, we have the $$1 \times N$$ curing rates vector $$\Delta=(\delta_1, \ldots, \delta_N)$$, with components $$\delta_z=\overline{\delta}_j$$ for all $$z \in V_j$$ and $$j=1, \ldots, n$$. In Section A.1 of the Appendix we shall discuss when and how it is possible to reduce the original system (2.2) to a system of $$n$$ differential equations, through the matrix $$Q^{\top}$$. 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 $$\pi=\left\{V_1, \ldots, V_n\right\}$$ be an equitable partition. Let $$A^{\top}$$ and $$Q^{\top}$$ be weighted matrices as in (4.2) and (4.4), respectively. Then, it holds: (i) $$(A^{\top}-D)S^{\top}=S^{\top}(Q^{\top}-\overline D)$$. (ii) For all $$\lambda \in \mathbb{C}$$ and all $$x \in \mathbb{C}^n$$ \begin{equation*} (Q^{\top}-\overline D)x= \lambda x \quad \text{if and only if} \quad (A^{\top}-D)S^{\top}x= \lambda S^{\top}x. \end{equation*} Now let us consider the system of $$N$$ differential equations (2.2). It is possible to extend [20, Theorem 4.1] to the case of directed graphs. Following that result, if we assume that at time $$t=0$$ 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 $$Q^{\top}$$. Hence, the reduced dynamical system writes \begin{align}\label{eqred} \frac{d\overline{p}_j(t)}{dt} &= (1-\overline{p}_j(t))\sum_{\substack{m=1 \\ m\neq j}}^n \!\! c^{\rm out}_{jm} \overline{p}_m(t) + c^{\rm out}_{jj}(1-\overline{p}_j(t))\overline{p}_j(t) - \overline{\delta_j} \overline{p}_j(t), \hskip15mm j=1,\ldots,n \end{align} (4.5) where $$\overline{p}_j(t)$$ is the infection probability of a node in the community $$j$$. We can prove that, in the case of a graph whose node set has an equitable partition, and regardless of initial conditions, the critical threshold for (2.2), applying Theorem 2.1, can be determined directly considering the reduced system (A.1). Proposition 4.3 The elements of the curing rates vector $$\Delta=( \delta_1, \,{\ldots}\,,\delta_{N})$$, that determine the critical threshold of (2.2), are identified by the elements of $$\overline \Delta=(\overline \delta_1, \,{\ldots}\, , \overline \delta_n)$$, in such a way that $$\delta_z=\overline \delta_j$$ for all $$z \in V_j$$, $$j=1, \ldots, n$$, for which \begin{equation}\label{qthr} r( Q^{\top} - \overline D) = 0, \end{equation} (4.6) where $$r$$ 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. 5. Optimization for networks with equitable partitions In this section, we consider a heterogeneous curing control per community. First, we assume that all nodes in community $$j$$ infect all nodes in community $$i$$ with the same infection rate, $$\beta_{V_i V_j}$$, and that $$\beta_{V_i V_j}={\beta}_{V_j V_i}$$, $$i,j=1, \ldots, n$$. In this case, the graph is undirected and the weights are symmetric, thus the quotient matrix $$Q$$ is symmetric and has real eigenvalues. Now let us consider the $$1 \times n$$ reduced curing rate vector $$\overline \Delta$$, the cost function writes \[ U(\overline \Delta)=\sum_{j=1}^n c_j k_j \overline{\delta}_j. \] Thus, $$U(\overline \Delta)$$ is the cost for curing all elements of each community $$j$$ at rate $$\overline{\delta}_j$$, where $$c_j >0$$, $$j=1, \ldots,n$$. We seek for the solution of the following Problem 5.1 (Eigenvalue constraint formulation) Find $$\overline \Delta \geq 0$$ which solves \begin{eqnarray*} && {\rm{minimize}} \qquad\qquad\quad U(\overline \Delta)\\ &&{\rm{subject\,\,to:}} \quad\;\; \lambda_1\!\big ( Q - {\rm diag}{(\overline \Delta)} \big ) \leq 0, \quad\;\; \overline \Delta \geq 0 \nonumber \end{eqnarray*} which also writes Problem 5.2 (SDP formulation) Find $$\overline \Delta \geq 0$$ which solves \begin{eqnarray*} && {\rm{minimize}} \qquad\qquad U(\overline \Delta)\nonumber \\ &&{\rm{subject\,\,to:}} \quad\;\; {\rm diag}{(\overline \Delta) - Q} \geq 0\nonumber \\ &&\quad\quad\quad\quad\quad\quad\;\; \overline \Delta \geq 0 \nonumber \end{eqnarray*} Theorem 3.3 guarantees the feasibility of the problem. The general case of equitable partitions introduced in Section 4.1, may not lead to a symmetric quotient matrix $$Q$$. However, we may consider suboptimal solutions—as explained in Section 3.1—obtained by applying the formulation of our optimization problem to the Hermitian part of $$Q$$. 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. 5.1 Two-level curing problem The state of the art for SDP solvers such as, for example the SDPT3 solver used for our numerical computation, provide solutions for moderate size graphs. Actually, the best known bound for the complexity of an $$\epsilon$$-solution attained with an interior point method is $$O( n^{3.5}\log(1/\epsilon))$$, where $$\epsilon$$ represents the accuracy [51]. 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 $${\varepsilon}$$-approximation of the optimal solution, with a complexity equal to $$O(\log(n)n^{3.3731} \log(1/\epsilon))$$ (see Theorem 5.7). Precisely, we consider only two possible levels of the nodes local curing rates, let us say $$\delta_0$$ and $$\delta_1$$, that are not fixed a priori. This situation fits well, for example in the case of a network where communities are of ‘two types’. Communities of the first type are eligible for curing rate $$\delta_0$$, whereas communities of the second type are eligible for curing rate $$\delta_1$$. For convenience, we define the former, central communities, and the latter, terminal communities. Such kind of configuration is suitable for a network that is, for example bipartite (where each node, e.g. represents a full-meshed community), or for an interconnected stars network, that is 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 [52], 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 [53]. 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, $$\pi_0=\left\{V_1^{0},\,{\ldots}\,,V_m^{0}\right\}$$ and $$\pi_1=\left\{V_1,\,{\ldots}\,,V_{m'}\right\}$$. We assume that the node set partition $$\pi=\pi_0 \cup \pi_1$$ is equitable. Let us introduce the curing matrix $$D={\rm diag}\left(\delta_0 \textbf{1}_{m}, \delta_1 \textbf{1}_{m^{'}}\right)$$ and define \[I^0_m =\begin{bmatrix} I_m& 0\\ 0 & 0 \end{bmatrix} , \qquad I^1_{m'} =\begin{bmatrix} 0 & 0\\ 0 & I_{m'} \end{bmatrix}\!, \] where $$I_m$$ is the identity matrix of order $$m$$. Then, we can write the SDP for the two-level curing rates, shortly the two-dimensional curing problem, as follows: Problem 5.3 (SDP two-dimensional formulation) Find $$\Delta_2=(\delta_0, \delta_1)$$ which solves \begin{eqnarray*} && {\rm{minimize}} \qquad\qquad\quad U(\Delta_2) \nonumber \\ &&{\rm{subject\,\,to:}} \quad\;\; {\delta_0 I^0_m + \delta_1 I^1_{m'} - Q} \geq 0\nonumber \\ &&\quad\quad\quad\quad\quad\quad\;\; \Delta_2 \geq 0 \nonumber \end{eqnarray*} The cost function is \[ U(\Delta_2)= \sum_{V_j \in \pi_0} k_j f_0(\delta_0) + \sum _{V_z \in \pi_1} k_z f_1(\delta_1), \] where $$f(\delta_0)=c_0 \delta_0$$, and $$f_1(\delta_1)= c_1 \delta_1$$, with $$ c_0,c_1 > 0$$, represent the effort to modify the recovery rate for nodes belonging to $$V_j \in \pi_0$$, and $$V_z \in \pi_1$$, respectively. In Section A.5 of the Appendix, we shall provide some simple examples for the optimal solution of the Problem (5.3). 5.2 Properties of the two-dimensional curing problem In the design of our algorithmic solution, we have leveraged on some basic properties of the two-dimensional curing problem. In order to do so, we need a few basic facts recalled next. Proposition 5.4 Let $$A$$ be an $$n\times n$$ symmetric, irreducible and non-negative matrix and let $$D={\rm diag}(\delta_1, \,{\ldots}\,,\delta_n)$$: i. if $$\delta_i=0$$ for some $$i=1,\ldots n$$, then $$\lambda_1(A-D) \geq 0$$; ii. The function $$(\delta_1, \,{\ldots}\,,\delta_n) \longmapsto$$$$\lambda_1(A-D)$$ is continuous; iii. $$\lambda_1(A-D)$$ is strictly decreasing in $$\delta_i$$, $$i=1,\ldots,n$$. Let us denote by $$\Gamma=\{(\delta_0,\delta_1)|\, \lambda_1({\rm diag}{(\delta_0\mathbf{1}_m,\delta_1 \mathbf{1}_{m'})} - Q)\geq 0\}$$ the feasibility region of Prob. 5.3: same argument of Theorem 3.3 let us state that it is non-empty; furthermore, the problem structure guarantees $$\Gamma$$ to be convex [31]. We indicate $$\Gamma_0$$ and $$\Gamma_1$$ the standard projections of $$\Gamma$$ onto the $$\delta_0$$-axis and the $$\delta_1$$-axis, respectively. Lemma 5.1 (Monotonicity) Let $$\phi:\delta_0\longmapsto \phi(\delta_0)$$ be the function that associates to $$\delta_0 \in \Gamma_0$$ the value $$\delta_1=\phi(\delta_0) \in \Gamma_1$$ such that $$\lambda_1( Q - {\rm diag}{(\delta_0\mathbf{1}_m,\delta_1 \mathbf{1}_{m'})} )= 0$$. Then, $$\phi$$ is decreasing. Proof. First, let us show that $$\phi$$ is a well-defined function over $$\Gamma_0$$. Let $$\delta_0\in \Gamma_0$$, because of feasibility, there exists $$\overline \delta_1$$, where $$(\delta_0,\overline \delta_1)\in\Gamma$$, such that $$\lambda_1( Q-{\rm diag}{(\delta_0\mathbf{1}_m, \overline \delta_1 \mathbf{1}_{m'})} )\leq 0$$. Furthermore, it holds $$\lambda_1(Q -{\rm diag}{(\delta_0\mathbf{1}_m,0\, \mathbf{1}_{m'})})\geq 0$$ by $$i)$$ of Proposition 5.4. Because of ii. and iii,in Proposition 5.4, we know that $$\lambda_1(Q-{\rm diag}{(\delta_0\mathbf{1}_m,\delta_1 \mathbf{1}_{m'})} )$$ is a continuous strictly decreasing function of $$\delta_1$$ over $$[0,\overline \delta_1]$$, so there exists one and only one value $$\delta_1 \in \Gamma_1$$ satisfying the definition of $$\phi$$. Let $$z>0$$ and assume that $$\phi(\delta_0+z)=\phi(\delta_0)+\zeta>\phi(\delta_0)$$, for some $$\zeta>0$$, that is that $$\phi$$ is not decreasing. From the definition of $$\phi$$ there exists $$0 \not = w \in \ker\Big(\!{\rm diag}{(((\delta_0+z)\mathbf{1}_m,\phi(\delta_0+z)\mathbf{1}_{m'})} - Q \Big )$$. Hence, we can write \begin{equation*} \begin{aligned}\label{eq:unique} &w^{\top} \big( Q - {\rm diag}\big(\delta_0\mathbf{1}_m,\phi(\delta_0)\mathbf{1}_{m'}\big)\big) w = w^{\top} {\rm diag}{(z\mathbf{1}_m,\zeta\mathbf{1}_{m'})} w\nonumber + w^{\top} \big(Q - {\rm diag}\big((\delta_0+z)\mathbf{1}_m,\phi(\delta_0+z)\mathbf{1}_{m'}\big) \big) w \nonumber \\ &\quad{} = w^{\top} {\rm diag}{((z\mathbf{1}_m,\zeta\mathbf{1}_{m'}))} w > 0, \end{aligned} \end{equation*} where the strict inequality holds because $${\rm diag}{(z\mathbf{1}_m,\zeta\mathbf{1}_{m'})} > 0$$. Since $$\lambda_1\Big( Q - {\rm diag}(\delta_0\mathbf{1}_m,\phi(\delta_0)\mathbf{1}_{m'} \Big )=0$$, this means that $$ Q - {\rm diag}(\delta_0\mathbf{1}_m,\phi(\delta_0)\mathbf{1}_{m'}) $$ 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 $$\Gamma$$. Theorem 5.5 (Compact search set) There exist two pairs $$(\delta_0^{\min},\delta_0^{\max})$$ and $$(\delta_1^{\min},\delta_1^{\max})$$ such that a solution $$\Delta_2^*=(\delta_0^*, \delta_1^*)$$ of Problem 5.3 belongs to a compact subset $$\Gamma' \subseteq [\delta_0^{\min}, \delta_0^{\max}] \times [\delta_1^{\min}, \delta_1^{\max}] $$. Proof. Let us define $$\hat{c}_0=\sum_{V_j \in \pi_0} c_0 k_j$$ and $$\hat{c}_1=\sum_{V_z \in \pi_1} c_1 k_z$$, then we write $$U_{l_{\max}}=\hat{c}_0 l_{max} + \hat{c}_1 l_{\max}$$, with $$l_{\max}$$ as in Theorem 3.3, and $$U^*=\hat c_0 \delta^*_0 + \hat c_1 \delta^*_1$$. Let us denote $$\Delta^{l_{\max}}_2=\left( l_{\max}, l_{\max}\right)$$, by Theorem 3.3, $$\Delta_2^{l_{\max}} \in \Gamma$$, hence $$U_{l_{\max}}\geq U^*$$ and, by defining set $$\Omega=\{(\delta_0, \delta_1):\hat c_0 \delta_0 + \hat{c_1}\delta_1\leq U_{l_{\max}})\}$$, it follows that $$(\delta_0^*,\delta_1^*)\in \Gamma'=\Gamma \cap \Omega$$; $$\Gamma'$$ is closed as intersection of closed sets. Now, feasibility conditions of Problem 5.3 require matrix $$ Q- \left(\delta_0 I^0_m + \delta_1 I^1_{m'}\right) $$ to be semidefinite negative. We define $$f(\delta_0)=\lambda_1\big( Q - \big(\delta_0 I^0_m + (\frac{U_{l_{\max}}-\hat c_0 \delta_0}{\hat c_1})I^1_{m'}\big) \big)$$: we have $$f(l_{\max})\leq 0$$ since $$( l_{\max}, l_{\max}) \in\Gamma$$ and $$f(0)\geq 0$$ by i. of Problem 5.4. By assertion ii. in Prop. 5.4, $$f(\delta_0)$$ is a continuous function. Hence, there exists $$\delta_0^{\min}$$ such that $$f(\delta_0^{\min})=0$$, and since $$\phi$$ is decreasing $$\phi(\delta_0^{\min})=\delta_1^{\max}$$. We can repeat the same reasoning by inverting the role of $$\delta_1$$ and $$\delta_0$$ defining $$g(\delta_1)=\lambda_1\big( Q - \big((\frac{U_{l_{\max}}-\hat c_1 \delta_1}{\hat c_0}) I^0_m +\delta_1 I^1_{m'}\big)\big)$$. Hence, we can assert that exists $$\delta_1^{\min}$$ such that $$g(\delta_1^{\min})=0$$ and $$\phi(\delta_1^{\min})=\delta_0^{\max}$$. Finally, by letting $$r: \hat c_0\delta_0 + \hat c_1 \delta_1 = U_{l_{\max}}$$, the points $$(\delta_0^{\min},\delta_1^{\max})$$ and $$(\delta_0^{\max},\delta_1^{\min})$$ belong to $$\partial \Gamma \cap r$$, that is they belong to $$\partial\Gamma'$$, so $$\Gamma' \subseteq [\delta_0^{\min}, \delta_0^{\max}] \times [\delta_1^{\min}, \delta_1^{\max}]$$, and consequently, being $$\Gamma'$$ closed, it is also compact. □ Remark 5.1 Theorem 5.5 allows us to identify an interval of the values of $$\delta_0$$ and $$\delta_1$$ where we can restrict the search of $$(\delta_0^*,\delta_1^*)$$. Since $$\Gamma' \subseteq [\delta_0^{\min}, \delta_0^{\max}] \times [\delta_1^{\min}, \delta_1^{\max}]$$ and $$(\delta_0^*, \delta_1^*) \in \Gamma'$$, then $$\delta_0^* \in [\delta_0^{\min}$$, $$\delta_0^{\max}]$$ and $$\delta_1^* \in [\delta_1^{\min}, \delta_1^{\max}]$$. This is one key property in the algorithmic search of the optimal solution presented in the following section. Finally, a direct proof that the optimal solution lies on $$\partial \Gamma'$$ follows: Corollary 5.1 A solution $$\Delta_2^*=(\delta_0^*, \delta_1^*)$$ of Problem 5.3 belongs to $$\partial \Gamma' \cap \Omega$$. Proof. Let us assume $$\Delta_2^*=(\delta_0^*,\delta_1^*) \in \Gamma' \setminus \partial \Gamma'$$. $$\Delta_2^*$$ is feasible, hence $$\lambda_1( Q - D)<0$$, with $$D={\rm diag}\left(\delta^*_0 \textbf{1}_{m}, \delta^*_1 \textbf{1}_{m'}\right)$$. From Prop. 5.4, again we can find $$0<\delta_1'<\delta_1^*$$ such that $$\lambda_1( Q - {\rm diag}(\delta_0^*\textbf{1}_{m},\delta_1'\textbf{1}_{m'}))=0$$, where, that is $$\Delta_2'=(\delta_0^*,\delta_1')\in \partial \Gamma'$$. But, $$U(\Delta_2^*)-U(\Delta_2')=\hat c_1(\delta_1^*-\delta_1')>0$$. Contradiction. □ 5.3 Bisection algorithm Table 1 reports on the pseudocode of algorithm OptimalThreshold2D: it solves the two-dimensional curing problem. It employs three additional functions LeftCorner (Table 2) RightCorner and, finally, BisectionThreshold (Table 3). Table 1 OptimalThreshold2D: solves the two-dimensional optimal curing problem via the bisection search Table 1 OptimalThreshold2D: solves the two-dimensional optimal curing problem via the bisection search Table 2 LeftCorner: identifies the left corner of $$\Gamma' \subseteq \Gamma$$ (Theorem 5.5); the pseudocode of the dual function $$(\delta_0^{\max},\delta_1^{\min}) =$$RightCorner$$(Q,c_0,c_1)$$ is omitted for the sake of space Table 2 LeftCorner: identifies the left corner of $$\Gamma' \subseteq \Gamma$$ (Theorem 5.5); the pseudocode of the dual function $$(\delta_0^{\max},\delta_1^{\min}) =$$RightCorner$$(Q,c_0,c_1)$$ is omitted for the sake of space Table 3 BisectionThreshold: given feasible $$\delta_0$$, finds $$\delta_1$$ such that $$(\delta_0,\delta_1)$$ lies on the frontier of the feasibility region Table 3 BisectionThreshold: given feasible $$\delta_0$$, finds $$\delta_1$$ such that $$(\delta_0,\delta_1)$$ lies on the frontier of the feasibility region LeftCorner identifies via bisection feasible point $$(\delta_0^{\min}, \delta_1^{\max})$$; the bisection search operated by LeftCorner—see proof of Theorem 5.5—is performed along values $$\delta_1=f(\delta_0)$$. The companion function RightCorner identifies the point $$(\delta_0^{\max}, \delta_1^{\min})$$; the pseudocode is omitted for the sake of space. Procedure isNegativeDefinite is the standard test for a real symmetric matrix $$A$$ to be negative definite; it requires to verify $${\mathrm{sgn}} \left(\det(A_k)\right)=(-1)^k$$ where $$A_k$$ is the $$k$$-th principal minor of $$A$$, that is the matrix obtained considering the first $$k$$ rows and columns only. Finally, the OptimalThreshold2D algorithm performs a bisection search based on a subgradient descent over the utility function $$U(\delta_0)=\hat c_0 \delta_0 + \hat c_1 \phi(\delta_0)$$. Remark 5.2 In Table 1, we have reported an implementation assuming the calculation of the subgradient $$\partial U$$ at each mid-point $$x$$. However, it is sufficient to evaluate the increment at a point $$x+\epsilon_1$$ within the feasibility region for some $$\epsilon_1>0$$: if $$U(x)<U(x+\epsilon_1)$$, then, due to convexity, the whole interval $$[x+\epsilon_1,+\infty)$$ can be discarded. Conversely, if $$U(x)>U(x+\epsilon_1)$$, then, due to convexity, the whole interval $$[0,x)$$ can be discarded during the search. This operation can be performed at a cost $$O(1)$$ when $$U(x)$$ and $$U(x+\epsilon_1)$$ are known, that is at the cost of two calls of BisectionThreshold. We note that REPEAT loop stops when $$\epsilon>\prod |\lambda_i|=|\det( Q - D)|>|\lambda_1|^n$$, that is when $$|\lambda_1|< (\epsilon)^{1/n}$$. Furthermore, the termination condition in BisectionThreshold, LeftCorner and RightCorner requires $$\Delta_2$$ to lie within the feasible region and the determinant to be smaller than $$\epsilon$$. Theorem 5.6 (Correctness) OptimalThreshold2D is an $$\epsilon$$-approximation of an optimal solution. Proof. The algorithm operates a bisection search for a global minimum of $$U(\Delta_2)=\hat c_1 \delta_0 + \hat c_1 \phi(\delta_0)$$, where $$U(\Delta_2)$$ is a convex function. Let $$V=U_{l_{\max}}$$ and $$\Delta_2^*$$ be the optimal solution: from the properties of the bisection search on (quasi-)convex functions [31, Chapter 4, pp. 145], the accuracy at step $$r= \lceil\log_2 ( V/\epsilon ) \rceil $$ of the algorithm is $$|U_r - U(\Delta_2^*)|< V\, 2^{-r} < \epsilon$$. □ Furthermore, we can characterize the computational complexity of the algorithm. Theorem 5.7 (Complexity) The time complexity of OptimalThreshold2D is $$O(n^{1+\ell} \log_2(n/\epsilon) )$$ where $$\ell=2.373$$. Proof. The number of iterations of the bisection search WHILE loop (lines $$1$$ to $$9$$ in Table 1) is $$O(\log_2(n/\epsilon))$$. This follows again from elementary properties of bisection search [31, Chapter 4, pp. 145]. In fact, the bisection search operates for $$0\leq U(\delta_0) \leq U_{l_{\max}}$$ and $$U_{l_{\max}}=l_{\max}(\hat{c}_0+\hat{c}_1)$$. Finally, indeed, $$l_{\max}\leq (n-1) \max_{i,j} q_{ij}$$. Same argument on the measure of the search intervals of BisectionThreshold, LeftCorner and RightCorner let us conclude that they require $$O(\log_2(n/\epsilon))$$ iterations of the REPEAT loop. Finally, test isNegativeDefinite appearing in Threshold2D, LeftCorner and Right-Corner requires the computation of $$n-1$$ determinants of the principal minors of $$A - D$$ at cost $$O(n^{1+\ell})$$. Here $$\ell$$ is the exponent for fast matrix multiplication [54]. In the case of the Coppersmith–Winograd algorithm for fast matrix multiplication, it holds $$\ell=2.373$$. □ 5.4 Numerical Results In this section we present the results of numerical experiments in the case of interconnected stars networks, a sample network is depicted in Fig. 2(b). In Fig. 3(a), we compare the ratio between the cost $$U_u$$ of the uniform curing rate vector, and the optimal cost $$U^*=U(\Delta^*)$$ obtained by solving the two-dimensional curing Problem 5.3, by means of the OptimalThreshold2D. The uniform curing rate vector is $$\Delta= \delta {\mathbf 1}_N$$, where $$\delta$$ is the value such that the threshold in (4.6) is attained. For this experiment, we consider that the infection spreads with rate $$\beta_{V_i^0 V_j^0}=\beta_0$$ among the central communities and with rate $$\beta_{V_i V_i^0}=\beta_1$$ between a central node and a node in its adjacent terminal community; moreover, we assume $$c_0=c_1=1$$. We consider that each terminal community has the same number of elements $$k$$. The computation is made for different values of $$k$$, for three different sample networks, with $$m=50$$ central nodes and $$m'=50$$ terminal communities. Sample networks differ for the average degree of the central nodes. Central nodes are connected as Erdős—Rényi random graphs with $$p=0.2$$, $$p=0.3$$, $$p=0.6$$, respectively. The plot confirms that a larger gain is obtained, in terms of costs, by two-dimensional curing policy versus a uniform approach, in particular, the larger the denser the network, namely, for larger $$p$$ in our samples. For the interconnected stars networks samples, in particular, we observe one order of magnitude gain in the cost function. We see that the advantage increases as the number of elements $$k$$ increases, with a $$\sqrt k$$ shaped ratio (A.9) as derived in closed form for the case $$m=m'=2$$. In Fig. 3(b), we have instead reported, only, on the optimal cost $$U^*$$ for different values of $$c_0$$ and $$c_1$$. In particular, we observe that the optimal cost appears to depend linearly on the community size $$k$$. Larger costs are incurred in the case when the coefficient $$c_0$$, related to the expenditure for the central nodes, is larger than $$c_1$$. This is in line with the fact that central communities are more connected than terminal communities, and consequently we need for more investments in such a way that the infection is kept subcritical. Fig. 3. View largeDownload slide (a) Ratio $$U_u/U^*$$ for increasing size $$k$$ of the terminal communities of interconnected star networks with $$m=m'=50$$. The three curves refer to networks where the central nodes are connected as Erdős—Rényi graphs generated for $$p=0.2$$, $$p=0.3$$ and $$p=0.6$$ respectively; $$\beta_0=1$$, $$\beta_1=0.3$$ and $$c_0=c_1=1$$. (b) Cost function $$U^*$$ for increasing dimension $$k$$ of the terminal communities, $$\beta_0=1$$, $$\beta_1=0.3$$. The curves refer to the case $$p=0.3$$ in the cases $$c_0=2c_1$$, $$c_0=c_1$$ and $$2 c_0=c_1$$, respectively. Fig. 3. View largeDownload slide (a) Ratio $$U_u/U^*$$ for increasing size $$k$$ of the terminal communities of interconnected star networks with $$m=m'=50$$. The three curves refer to networks where the central nodes are connected as Erdős—Rényi graphs generated for $$p=0.2$$, $$p=0.3$$ and $$p=0.6$$ respectively; $$\beta_0=1$$, $$\beta_1=0.3$$ and $$c_0=c_1=1$$. (b) Cost function $$U^*$$ for increasing dimension $$k$$ of the terminal communities, $$\beta_0=1$$, $$\beta_1=0.3$$. The curves refer to the case $$p=0.3$$ in the cases $$c_0=2c_1$$, $$c_0=c_1$$ and $$2 c_0=c_1$$, respectively. In Table 4, we compare the performance of OptimalThreshold2D with an SDP solver, namely the SDPT3 solver [34]. The SDPT3 solver generates a solution using a primal–dual interior-point algorithm which leverages on the infeasible path-following paradigm. As reported in Table 4, when the solver is applied to Problem 5.3, we denote the corresponding solution as SDPT3(2D). For the sake of comparison, we have reported also on the optimal solution derived with the same solver when curing rates are optimized per node (Problem 3.2), and we refer to this solution as SDPT3. The solution is provided on a graph with $$m=m'=50$$ and $$c_0=c_1=1$$, for increasing values of the terminal community dimension $$k$$. We can observe that for the interconnected stars network, in the case of two infection rate levels, SDPT3, SDPT3(2D) and OptimalThreshold2D provide similar values. This result suggests that, in this particular case, there is no advantage to treat each node with different curing policies: the general solution obtained using SDPT3 has similar performance as the one generated solving the two-dimensional formulation of the problem, that is by using OptimalThreshold2D or SDPT3(2D) in the two-dimensional case. We observe also that the curing rate of terminal nodes appears insensitive to the increase of the terminal communities size $$k$$, as it can be seen, with a direct computation, in the Section A.5 of the Appendix for the case with $$m=1$$ and $$m=2$$, respectively. In Table 5, same performance evaluation has been reported for the same sample graph and the same cost coefficients, but studying the case with more than two infection rate levels; specifically a central node can eventually infect each of its adjacent central nodes with a different infection rate, also the infection rate between a central node and a terminal community can vary from a subgraph to another. The infection rates are generated as uniform random variables in the interval $$(0.1,1.9)$$ and $$(0.03,0.57)$$ for the speed of infection between central communities and between a central node and a terminal community, respectively. Table 4 Performance of OptimalThreshold2D, SDPT3 (2D) and SDPT3. The graph considered is an interconnected star network with $$m=m'=50$$, where the connection between the central nodes are represented by a Erdős—Rényi graph generated with $$p=0.2$$; $$c_0=c_1=1$$ and the values of the weights are $$\beta_0=1$$ and $$\beta_1=0.3$$. In the SDPT3 case, the values of $$\delta_i^*$$, $$i=0,1$$, represent the averaged value of the node-specific curing rates over a community OptimalThreshold2D SDPT3 (2D) SDPT3 k $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$ \delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ 10 0.8057 13.116 0.29979 0.80572 13.1144 0.3 0.772 12.44 0.3 20 1.1057 16.1163 0.29984 1.1057 16.1144 0.3 1.072 15.44 0.3 30 1.4056 19.1222 0.29965 1.4057 19.1145 0.3 1.372 18.44 0.3 40 1.7055 22.1294 0.29951 1.7057 22.1144 0.3 1.672 21.44 0.3 50 2.0053 25.1354 0.29943 2.0057 25.1144 0.3 1.972 24.44 0.3 60 2.3052 28.1414 0.29936 2.3057 28.1144 0.3 2.272 27.44 0.3 70 2.6049 31.1578 0.29916 2.6057 31.1145 0.3 2.572 30.44 0.3 80 2.9047 34.1792 0.29893 2.9057 34.1144 0.3 2.872 33.4401 0.3 90 3.2043 37.1929 0.29882 3.2057 37.1144 0.3 3.172 36.4402 0.3 100 3.5039 40.1319 0.29947 3.5057 40.1144 0.3 3.472 39.44 0.3 OptimalThreshold2D SDPT3 (2D) SDPT3 k $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$ \delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ 10 0.8057 13.116 0.29979 0.80572 13.1144 0.3 0.772 12.44 0.3 20 1.1057 16.1163 0.29984 1.1057 16.1144 0.3 1.072 15.44 0.3 30 1.4056 19.1222 0.29965 1.4057 19.1145 0.3 1.372 18.44 0.3 40 1.7055 22.1294 0.29951 1.7057 22.1144 0.3 1.672 21.44 0.3 50 2.0053 25.1354 0.29943 2.0057 25.1144 0.3 1.972 24.44 0.3 60 2.3052 28.1414 0.29936 2.3057 28.1144 0.3 2.272 27.44 0.3 70 2.6049 31.1578 0.29916 2.6057 31.1145 0.3 2.572 30.44 0.3 80 2.9047 34.1792 0.29893 2.9057 34.1144 0.3 2.872 33.4401 0.3 90 3.2043 37.1929 0.29882 3.2057 37.1144 0.3 3.172 36.4402 0.3 100 3.5039 40.1319 0.29947 3.5057 40.1144 0.3 3.472 39.44 0.3 Table 4 Performance of OptimalThreshold2D, SDPT3 (2D) and SDPT3. The graph considered is an interconnected star network with $$m=m'=50$$, where the connection between the central nodes are represented by a Erdős—Rényi graph generated with $$p=0.2$$; $$c_0=c_1=1$$ and the values of the weights are $$\beta_0=1$$ and $$\beta_1=0.3$$. In the SDPT3 case, the values of $$\delta_i^*$$, $$i=0,1$$, represent the averaged value of the node-specific curing rates over a community OptimalThreshold2D SDPT3 (2D) SDPT3 k $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$ \delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ 10 0.8057 13.116 0.29979 0.80572 13.1144 0.3 0.772 12.44 0.3 20 1.1057 16.1163 0.29984 1.1057 16.1144 0.3 1.072 15.44 0.3 30 1.4056 19.1222 0.29965 1.4057 19.1145 0.3 1.372 18.44 0.3 40 1.7055 22.1294 0.29951 1.7057 22.1144 0.3 1.672 21.44 0.3 50 2.0053 25.1354 0.29943 2.0057 25.1144 0.3 1.972 24.44 0.3 60 2.3052 28.1414 0.29936 2.3057 28.1144 0.3 2.272 27.44 0.3 70 2.6049 31.1578 0.29916 2.6057 31.1145 0.3 2.572 30.44 0.3 80 2.9047 34.1792 0.29893 2.9057 34.1144 0.3 2.872 33.4401 0.3 90 3.2043 37.1929 0.29882 3.2057 37.1144 0.3 3.172 36.4402 0.3 100 3.5039 40.1319 0.29947 3.5057 40.1144 0.3 3.472 39.44 0.3 OptimalThreshold2D SDPT3 (2D) SDPT3 k $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$ \delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ 10 0.8057 13.116 0.29979 0.80572 13.1144 0.3 0.772 12.44 0.3 20 1.1057 16.1163 0.29984 1.1057 16.1144 0.3 1.072 15.44 0.3 30 1.4056 19.1222 0.29965 1.4057 19.1145 0.3 1.372 18.44 0.3 40 1.7055 22.1294 0.29951 1.7057 22.1144 0.3 1.672 21.44 0.3 50 2.0053 25.1354 0.29943 2.0057 25.1144 0.3 1.972 24.44 0.3 60 2.3052 28.1414 0.29936 2.3057 28.1144 0.3 2.272 27.44 0.3 70 2.6049 31.1578 0.29916 2.6057 31.1145 0.3 2.572 30.44 0.3 80 2.9047 34.1792 0.29893 2.9057 34.1144 0.3 2.872 33.4401 0.3 90 3.2043 37.1929 0.29882 3.2057 37.1144 0.3 3.172 36.4402 0.3 100 3.5039 40.1319 0.29947 3.5057 40.1144 0.3 3.472 39.44 0.3 Table 5 Performance of OptimalThreshold2D, SDPT3 (2D) and SDPT3, sample graphs are obtained from the same graph used in Table 4 and $$c_0=c_1=1$$. The infection rates are generated as uniform random variables in the interval $$(0.1,1.9)$$ and $$(0.03,0.57)$$, for rates between central communities and between a central node and a terminal community respectively. In the SDPT3 case, the values of $$\delta_i^*$$, $$i=0,1$$, represent the averaged value of the node-specific curing rates over a community OptimalThreshold2D SDPT3 (2D) SDPT3 k $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$ \delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ 10 1.9963 34.1309 0.57945 1.9963 34.1309 0.57945 1.9379 33.4102 0.53481 20 2.5603 39.3993 0.5903 2.5603 39.3993 0.5903 2.4307 38.3384 0.51382 30 3.1665 44.9001 0.61431 3.1666 44.9001 0.61431 2.9536 43.5669 0.51682 40 3.7954 50.6431 0.63164 3.7957 50.6431 0.63164 3.4781 48.8125 0.51876 50 4.4332 56.4137 0.64499 4.4336 56.4137 0.64499 4.0279 54.3098 0.52495 60 5.1272 62.627 0.66528 5.1278 62.627 0.66528 4.6049 60.0804 0.53364 70 5.8175 69.0209 0.67612 5.8184 69.0209 0.67612 5.1595 65.6263 0.53663 80 6.4779 74.8633 0.68369 6.4792 74.8633 0.68369 5.6708 70.7388 0.53346 90 7.1277 80.4688 0.68984 7.1298 80.4688 0.68984 6.1686 75.717 0.5295 100 7.8361 87.1312 0.6959 7.839 87.1312 0.6959 6.7182 81.2132 0.53151 OptimalThreshold2D SDPT3 (2D) SDPT3 k $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$ \delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ 10 1.9963 34.1309 0.57945 1.9963 34.1309 0.57945 1.9379 33.4102 0.53481 20 2.5603 39.3993 0.5903 2.5603 39.3993 0.5903 2.4307 38.3384 0.51382 30 3.1665 44.9001 0.61431 3.1666 44.9001 0.61431 2.9536 43.5669 0.51682 40 3.7954 50.6431 0.63164 3.7957 50.6431 0.63164 3.4781 48.8125 0.51876 50 4.4332 56.4137 0.64499 4.4336 56.4137 0.64499 4.0279 54.3098 0.52495 60 5.1272 62.627 0.66528 5.1278 62.627 0.66528 4.6049 60.0804 0.53364 70 5.8175 69.0209 0.67612 5.8184 69.0209 0.67612 5.1595 65.6263 0.53663 80 6.4779 74.8633 0.68369 6.4792 74.8633 0.68369 5.6708 70.7388 0.53346 90 7.1277 80.4688 0.68984 7.1298 80.4688 0.68984 6.1686 75.717 0.5295 100 7.8361 87.1312 0.6959 7.839 87.1312 0.6959 6.7182 81.2132 0.53151 Table 5 Performance of OptimalThreshold2D, SDPT3 (2D) and SDPT3, sample graphs are obtained from the same graph used in Table 4 and $$c_0=c_1=1$$. The infection rates are generated as uniform random variables in the interval $$(0.1,1.9)$$ and $$(0.03,0.57)$$, for rates between central communities and between a central node and a terminal community respectively. In the SDPT3 case, the values of $$\delta_i^*$$, $$i=0,1$$, represent the averaged value of the node-specific curing rates over a community OptimalThreshold2D SDPT3 (2D) SDPT3 k $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$ \delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ 10 1.9963 34.1309 0.57945 1.9963 34.1309 0.57945 1.9379 33.4102 0.53481 20 2.5603 39.3993 0.5903 2.5603 39.3993 0.5903 2.4307 38.3384 0.51382 30 3.1665 44.9001 0.61431 3.1666 44.9001 0.61431 2.9536 43.5669 0.51682 40 3.7954 50.6431 0.63164 3.7957 50.6431 0.63164 3.4781 48.8125 0.51876 50 4.4332 56.4137 0.64499 4.4336 56.4137 0.64499 4.0279 54.3098 0.52495 60 5.1272 62.627 0.66528 5.1278 62.627 0.66528 4.6049 60.0804 0.53364 70 5.8175 69.0209 0.67612 5.8184 69.0209 0.67612 5.1595 65.6263 0.53663 80 6.4779 74.8633 0.68369 6.4792 74.8633 0.68369 5.6708 70.7388 0.53346 90 7.1277 80.4688 0.68984 7.1298 80.4688 0.68984 6.1686 75.717 0.5295 100 7.8361 87.1312 0.6959 7.839 87.1312 0.6959 6.7182 81.2132 0.53151 OptimalThreshold2D SDPT3 (2D) SDPT3 k $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$ \delta_0^*$$ $$k \cdot \delta_1^*$$ $$U^*(10^3)$$ $$\delta_0^*$$ $$k \cdot \delta_1^*$$ 10 1.9963 34.1309 0.57945 1.9963 34.1309 0.57945 1.9379 33.4102 0.53481 20 2.5603 39.3993 0.5903 2.5603 39.3993 0.5903 2.4307 38.3384 0.51382 30 3.1665 44.9001 0.61431 3.1666 44.9001 0.61431 2.9536 43.5669 0.51682 40 3.7954 50.6431 0.63164 3.7957 50.6431 0.63164 3.4781 48.8125 0.51876 50 4.4332 56.4137 0.64499 4.4336 56.4137 0.64499 4.0279 54.3098 0.52495 60 5.1272 62.627 0.66528 5.1278 62.627 0.66528 4.6049 60.0804 0.53364 70 5.8175 69.0209 0.67612 5.8184 69.0209 0.67612 5.1595 65.6263 0.53663 80 6.4779 74.8633 0.68369 6.4792 74.8633 0.68369 5.6708 70.7388 0.53346 90 7.1277 80.4688 0.68984 7.1298 80.4688 0.68984 6.1686 75.717 0.5295 100 7.8361 87.1312 0.6959 7.839 87.1312 0.6959 6.7182 81.2132 0.53151 As seen there, by curing nodes with different curing policies it is possible to attain lower costs at larger values of $$k$$. This effect is depicted also in Fig. 4. In particular, again, the two-dimensional curing rates output of SDPT3 (2D) and OptimalThreshold2D show similar performance and the optimal curing rate of terminal communities appears insensitive to the increase of the terminal communities size. We observe that the relative advantage of the SDPT3 tends to increase with the size of the terminal communities. Fig. 4. View largeDownload slide Details of the costs in the case of uniform distribution of infection rates (see Table 5). Fig. 4. View largeDownload slide Details of the costs in the case of uniform distribution of infection rates (see Table 5). In the final set of experiments we study the case of a complete bipartite graph. We consider that the community whose curing rate is $$\delta_0$$, to which we refer to as the central community, has a fixed dimension $$k_0=50$$; instead, for the so-called terminal community, with $$\delta_1$$ curing rate, we consider increasing size $$k_1=1, 50, 100, 150, 200$$. In Fig. 5(a), we report on the ratio between the cost obtained by using the uniform curing policy, namely $$U_u$$, and the optimal cost $$U^*$$ obtained by means of the OptimalThreshold2D, in the case of equal coefficients $$c_0=c_1=1$$. As expected, we can see that when $$k_1=1$$ we obtain an advantage in the use of the two-level curing strategy. Clearly, when the two communities have the same size there is no difference between the two costs; the ratio starts to grow again as the asymmetry in terms of communities dimensions, starts to increase again. Fig. 5. View largeDownload slide (a) Ratio $$U_u/U^*$$ in the case of complete bipartite graphs for increasing size $$k_1=1, 50, 100, 150, 200$$ of the terminal community, and fixed size $$k_0=50$$ of the central community: $$\beta=1$$, $$c_0=c_1=1$$. (b) Cost function $$U^*$$ for increasing dimension $$k_1$$ of the terminal community and fixed size $$k_0=50$$, $$\beta=1$$. The two curves refer to the cases $$c_0=c_1=1$$ and $$c_1=4c_0$$, respectively, in such a way that $$c_0+c_1=2$$. Fig. 5. View largeDownload slide (a) Ratio $$U_u/U^*$$ in the case of complete bipartite graphs for increasing size $$k_1=1, 50, 100, 150, 200$$ of the terminal community, and fixed size $$k_0=50$$ of the central community: $$\beta=1$$, $$c_0=c_1=1$$. (b) Cost function $$U^*$$ for increasing dimension $$k_1$$ of the terminal community and fixed size $$k_0=50$$, $$\beta=1$$. The two curves refer to the cases $$c_0=c_1=1$$ and $$c_1=4c_0$$, respectively, in such a way that $$c_0+c_1=2$$. In Fig. 5(b), we compare the effect of having two different cost coefficients, precisely $$c_1=4 c_0$$, with the case where they are the same, namely $$c_0=c_1=1$$, considering that their total amount is given and fixed a priori, that is $$c_0+c_1=2$$. A bias in the cost coefficient towards terminal community has a beneficial consequence because the optimal cost is lower than the case of equal coefficients, the advantage increases the larger the size of the terminal community. However, it is interesting to note that we would have obtained the same optimal cost if we interchanged the two coefficients, namely $$c_0=4c_1$$ (see Section A.5). Basically, due to the specific topology of the network, we have an advantage in considering different cost coefficients for the two communities, instead of having them equal, if their total sum is fixed. 6. Conclusions We have studied the problem of finding a non-uniform allocation of curing resources, within a networked population, at the minimum cost possible to avoid the epidemic from persisting indefinitely in the network. We have considered a mean-field approximation of an SIS model to study the diffusion of epidemics over a directed weighted graphs, capturing the possible asymmetric interactions, and the heterogeneity in the contagiousness. We have reported on the necessary and sufficient conditions for the extinction of the epidemics. These conditions are related to the sign of the stability modulus of a matrix encoding for the network structure and for the parameters of the model. Thus, such stability modulus represents the epidemic threshold of our system. Consequently, we have formulated a convex optimization problem for determining a cost-optimal curing policy, via a SDP approach, involving, the spectral properties of the network. Then, we have specialized the theory to the case of equitable partitions, in order to model heterogeneous community networks that possess a certain degree of regularity in their connectivity; this choice has been motivated since communities are relevant non-trivial topological feature of complex networks, that often have a certain regularity in their structure. Thus, in this case, we were able to reduce the dimensionality of our optimization problem, that is useful since the size of many real-networks poses limitations in investigating their spectral properties. At last, we have discussed on the special case of a two-dimensional curing policy, that can reflect, for example the case of different policy decisions for two different kinds of individual units (male and female, younger and elders, small villages and cities, firewall/gateways or clients in an enterprise network, etc.,$${\ldots}$$). With respect to this problem we have proposed an $$\epsilon$$-approximation algorithm with polynomial complexity in the input size. Funding European Commission within the framework of the CONGAS project (FP7-ICT-2011-8-317672) in part. Footnotes 1 As suggested in [2], see [32] for a method to check if a matrix is symmetrizable, and, in case, the way to choose the diagonal matrix to achieve symmetry. References 1. Pastor-Satorras R. , Castellano C. , Van Mieghem P. & Vespignani A. ( 2015 ) Epidemic processes in complex networks. Rev. Modern Phys. , 87 , 925 . Google Scholar Crossref Search ADS 2. Wan Y. , Roy S. & Saberi A. ( 2008 ) Designing spatially heterogeneous strategies for control of virus spread. Syst. Biol., IET , 2 , 184 – 201 . Google Scholar Crossref Search ADS 3. Prakash B. A. , Adamic L. , Iwashyna T. , Tong H. & Faloutsos C. ( 2013 ) Fractional immunization in networks. In Proceedings of the 2013 SIAM International Conference on Data Mining (pp. 659 – 667 ). Austin, TX, USA : Society for Industrial and Applied Mathematics . 4. Borgs C. , Chayes J. , Ganesh A. & Saberi A. ( 2010 ) How to distribute antidote to control epidemics. Random Struct. Algorithms , 37 , 204 – 222 . 5. Gourdin E. , Omic J. & Van Mieghem P. ( 2011 ) Optimization of network protection against virus spread. DRCN . ( Jajszczyk A. Cholda P. & Helvik B. E. ed.), Piscataway, NJ, USA : IEEE Society , pp. 86 – 93 . 6. Sahneh F. D. & Scoglio C. M. ( 2012 ) Optimal information dissemination in epidemic networks. In IEEE 51st Annual Conference on Decision and Control (CDC), 2012 . IEEE , pp. 1657 – 1662 . 7. Van Mieghem P. , Omic J. & Kooij R. ( 2009 ) Virus spread in networks. IEEE/ACM Trans. Netw. 17 , 1 – 14 . Google Scholar Crossref Search ADS 8. Preciado V. M. , Zargham M. , Enyioha C. , Jadbabaie A. & Pappas G. ( 2013 ) Optimal vaccine allocation to control epidemic outbreaks in arbitrary networks. In Decision and Control (CDC) , 2013 IEEE 52nd Annual Conference on IEEE , pp. 7486 – 7491 . 9. Preciado V. M. , Zargham M. , Enyioha C. , Jadbabaie A. & Pappas G. ( 2014 ) Optimal resource allocation for network protection against spreading processes. IEEE Trans. Control Netw. Syst. , 1 , 99 – 108 . Google Scholar Crossref Search ADS 10. Enyioha C. , Jadbabaie A. , Preciado V. M. & Pappas G. ( 2015 ) Distributed resource allocation for epidemic control. European Control conference , arXiv preprint arXiv:1501.01701 . 11. Drakopoulos K. , Ozdaglar A. & Tsitsiklis J. N. ( 2016 ) When is a network epidemic hard to eliminate? Math. Operations Research , 42 , 1 – 14 . Google Scholar Crossref Search ADS 12. Drakopoulos K. , Ozdaglar A. & Tsitsiklis J. N. ( 2015 ) A lower bound on the performance of dynamic curing policies for epidemics on graphs. In Decision and Control (CDC) , 2015 IEEE 54th Annual Conference on (pp. 3560 – 3567 ). IEEE . 13. Boccaletti S. , Latora V. , Moreno Y. , Chavez M. & Hwang D.-U. ( 2006 ) Complex networks: structure and dynamics. Phys. Rep. , 424 , 175 – 308 . Google Scholar Crossref Search ADS 14. Schwenk A. J. ( 1974 ) Computing the characteristic polynomial of a graph. In Graphs and Combinatorics ( Bari R. A. & Harary F. eds), Berlin, Heidelberg : Springer , pp. 153 – 172 . 15. Godsil C. ( 1980 ) Feasibility conditions for the existence of walk-regular graphs. Linear Algebra Appl. , 30 , 15 – 61 . Google Scholar Crossref Search ADS 16. Mugnolo D. ( 2014 ) Semigroup Methods for Evolution Equations on Networks . Cham, Switzerland : Springer . 17. Watts D. J. , Muhamad R. , Medina D. C. , & Dodds P. S. ( 2005 ) Multiscale, resurgent epidemics in a hierarchical metapopulation model. Proc. Natl. Acad. Sci. USA , 102 , 11157 – 11162 . Google Scholar Crossref Search ADS 18. May R. M. & Anderson R. M. ( 1984 ) Spatial heterogeneity and the design of immunization programs. Math. Biosci. , 72 , 83 – 111 . Google Scholar Crossref Search ADS 19. Riley S. , Fraser C. , Donnelly C. A. , Ghani A. C. , Abu-Raddad L. J. , Hedley A. J. , Leung G. M. , Ho L. M. , Lam T. H. , Thach T. Q. & Chau P. ( 2003 ) Transmission dynamics of the etiological agent of sars in hong kong: impact of public health interventions. Science , 300 , 1961 – 1966 . Google Scholar Crossref Search ADS PubMed 20. Bonaccorsi S. , Ottaviano S. , Mugnolo D. & De Pellegrini. F. ( 2015 ) Epidemic outbreaks in networks with equitable or almost-equitable partitions. SIAM J. Appl. Math. , 75 , 2421 – 2443 . Google Scholar Crossref Search ADS 21. Sahneh F. D. , Scoglio C. , & Van Mieghem P. ( 2013 ) Generalized epidemic mean-field model for spreading processes over multilayer complex networks. IEEE/ACM Trans. Netw. , 21 , 1609 – 1620 . Google Scholar Crossref Search ADS 22. Van Mieghem P. ( 2013 ) Decay towards the overall-healthy state in sis epidemics on networks. arXiv preprint arXiv:1310.3980 . 23. Draief M. & Massouli L. ( 2010 ) Epidemics and Rumours in Complex Networks . New York : Cambridge University Press . 24. Van Mieghem P. ( 2016 ) Approximate formula and bounds for the time-varying susceptible-infected-susceptible prevalence in networks. Phys. Rev. E , 93 , 052312 . Google Scholar Crossref Search ADS PubMed 25. Van Mieghem P. , Sahneh F. D. & Scoglio C. ( 2014 ) An upper bound for the epidemic threshold in exact Markovian SIR and SIS epidemics on networks. In Decision and Control (CDC) , 2014 IEEE 53rd Annual Conference on (pp. 6228 – 6233 ). IEEE . 26. Cator E. & Van Mieghem P. ( 2014 ) Nodal infection in Markovian SIS and SIR epidemics on networks are non-negatively correlated. Phys. Rev. E , 89 , 052802 . Google Scholar Crossref Search ADS 27. Van Mieghem P. & Omic J. ( 2013 ) In-homogeneous virus spread in networks. arxiv: 1306.2588 . 28. Lajmanovich A. & Yorke J. A. ( 1976 ) A deterministic model for Gonorrhea in a non-homogeneous population. Math. Biosci. , 28 , 221 – 236 . Google Scholar Crossref Search ADS 29. Hupert N. , Cuomo J. , Callahan M. A. , Mushlin A. I. & Morse S. S. ( 2004 ) Community-based mass prophylaxis: a planning guide for public health preparedness . Diane Publishing Company. US Department of Health and Human Services, Agency for Healthcare Research and Quality . 30. Boyd S. P. ( 1994 ) Semidefinite programming. SIAM Rev. , 38 , 49 – 95 . 31. Boyd S. P. & Vandenberghe L. ( 2004 ) Convex Optimization . New York, USA : Cambridge University Press . 32. Berman A. & Plemmons R. J. ( 1994 ) Nonnegative Matrices in the Mathematical Sciences . SIAM . 33. Zhang F. ( 2011 ) Matrix Theory: Basic Results and Techniques . New York : Springer Science & Business Media . 34. Tütüncü R. H. , Toh K. C. & Todd M. J. ( 2003 ) Solving semidefinite-quadratic-linear programs using SDPT3. Math. Program. , 95 , 189 – 217 . Google Scholar Crossref Search ADS 35. Ball F. , Britton T. , House T. , Isham V. , Mollison D. , Pellis L. & Scalia Tomba G. ( 2015 ) Seven challenges for metapopulation models of epidemics, including households models. Epidemics , 10 , 63 – 67 . Google Scholar Crossref Search ADS PubMed 36. Hanski I. & Ovaskainen O. ( 2003 ) Metapopulation theory for fragmented landscapes. Theoret. Popul. Biol. , 64 , 119 – 127 . Google Scholar Crossref Search ADS 37. Masuda N. ( 2010 ) Effects of diffusion rates on epidemic spreads in metapopulation networks. N. J. Phys. , 12 , 093009 . Google Scholar Crossref Search ADS 38. Ball F. , Mollison D. & Scalia-Tomba G. ( 1997 ) Epidemics with two levels of mixing. Ann. Appl. Probab. , 7 , 46 – 89 . Google Scholar Crossref Search ADS 39. Ross J. V. , House T. & Keeling M. J. ( 2010 ) Calculation of disease dynamics in a population of households. PLoS One , 5 , e9666, 2010 . Google Scholar Crossref Search ADS 40. Ball F. & Neal P. ( 2002 ) A general model for stochastic sir epidemics with two levels of mixing. Math. Biosci. , 180 , 73 – 102 . Google Scholar Crossref Search ADS PubMed 41. Pellis L. , Ball F. & Trapman P. ( 2012 ) Reproduction numbers for epidemic models with households and other social structures. I. definition and calculation of r 0. Math. Biosci. , 235 , 85 – 97 . Google Scholar Crossref Search ADS PubMed 42. Ball F. , Sirl D. & Trapman P. ( 2008 ) Network epidemic models with two levels of mixing. Math. Biosci. , 212 , 69 – 87 . Google Scholar Crossref Search ADS PubMed 43. Frank B. , David S. & Pieter T. et al. ( 2009 ) Threshold behaviour and final outcome of an epidemic on a random network with household structure. Adv. Appl. Probab. , 41 , 765 – 796 . Google Scholar Crossref Search ADS 44. Wang H. , Li Q. , D’Agostino G. , Havlin S. , Eugene Stanley H. & Van Mieghem P. ( 2013 ) Effect of the interconnected network structure on the epidemic threshold. Phys. Rev. E , 88 , 022801 . Google Scholar Crossref Search ADS 45. Bonaccorsi S. , Ottaviano S. , De Pellegrini F. , Socievole A. & Van Mieghem P. ( 2014 ) Epidemic outbreaks in two-scale community networks. Phys. Rev. E , 90 , 012810 . Google Scholar Crossref Search ADS 46. Girvan M. & Newman M. E. J. ( 2002 ) Community structure in social and biological networks. Proc. Natl. Acad. Sci. USA , 99 , 7821 – 7826 . Google Scholar Crossref Search ADS 47. Stewart I. , Golubitsky M. & Pivato M. ( 2003 ) Symmetry groupoids and patterns of synchrony in coupled cell networks. SIAM J. Appl. Dyn. Syst. , 2 , 609 – 646 . Google Scholar Crossref Search ADS 48. Golubitsky M. , Stewart I. & Török A. ( 2005 ) Patterns of synchrony in coupled cell networks with multiple arrows. SIAM J. Appl. Dyn. Syst. , 4 , 78 – 100 . Google Scholar Crossref Search ADS 49. Rahmani A. , Ji M. , Mesbahi M. & Egerstedt M. ( 2009 ) Controllability of multi-agent systems from a graph-theoretic perspective. SIAM J. Control Optim. , 48 , 162 – 186 . Google Scholar Crossref Search ADS 50. Aguilar C. O. & Gharesifard B. ( 2016 ) On almost equitable partitions and network controllability. In American Control Conference (ACC), 2016 . IEEE, pp. 179 – 184 . 51. Nemirovski A. ( 2007 ) Advances in convex optimization: conic programming. In Proceedings of International Congress of Mathematicians. pp. 413 – 444 . 52. Barabási A.-L. & Albert R. ( 1999 ) Emergence of scaling in random networks. Science , 286 , 509 – 512 . Google Scholar Crossref Search ADS PubMed 53. Li C. , van de Bovenkamp R. & Van Mieghem P. ( 2012 ) Susceptible-infected-susceptible model: a comparison of n-intertwined and heterogeneous mean-field approximations. Phys. Rev. E , 86 , 026116 . Google Scholar Crossref Search ADS 54. Aho A. V. , Hopcroft J. E. & Ullman J. D. ( 1974 ) The Design and Analysis of Computer Algorithms . Boston, USA : Addison-Wesley . 55. Horn R. A. & Johnson C. R. ( 2012 ) Matrix Analysis . New York, NY, USA : Cambridge University Press . 56. Varga R. S. ( 2009 ) Matrix Iterative Analysis , vol. 27. Springer Science & Business Media . Appendix A.1 Dimensionality reduction of the dynamical system (2.2) Let us define $$q^{\top}_{jm}$$ as the element of $$Q^{\top}$$ in position $$(i,j)$$. We know that \[ q^{\top}_{jm}= \frac{k_j c^{\rm out}_{jm}}{\sqrt{k_j k_m}}, \] hence \[ c^{\rm out}_{jm}= \frac{\sqrt{k_j k_m}}{k_j} q^{\top}_{jm}=\left(\frac{k_m}{k_j}\right)^{1/2} q^{\top}_{jm}. \] Thus we can write (4.5) in the following matrix form \begin{equation}\label{hetQ12} \frac{d \overline{P}(t)}{dt}=\left( \tilde{Q}-\overline{D}\right)\overline{P}(t)-\mathop{\rm diag}(\overline{P}(t))\tilde{Q}\overline{P}(t), \end{equation} (A.1) where $$\widetilde Q= \operatorname{diag}\left(\frac{1}{\sqrt{k_j}}\right) Q^{\top}\operatorname{diag} (\sqrt{k_j})$$. It is immediate that $$\sigma(Q^{\top})= \sigma(\tilde{Q})$$. By [20, Corollary 4.2], irrespective of the initial conditions of nodes in the same community, it is sufficient to compute the positive steady-state vector $$\overline{P}_{\infty}$$ of the reduced system (A.1) to obtain that of the original system (2.2). Indeed, as time elapses all nodes in the same community tend to have the same infection probabilities, thus the components of the steady-state vector $$P_{\infty}$$ corresponding to nodes in the same community are equal. A.2 Proof of Lemma 4.1 Proof. i. We first prove that $$A^{\top}S^{\top}=S^{\top}Q^{\top}$$. In fact, if $$i \in V_h$$, it holds \begin{eqnarray} (A^{\top}S^{\top})_{i,j}&=& \frac{c^{\rm out}_{hj}}{\sqrt{k_j}}, \end{eqnarray} (A.2) \begin{eqnarray} (S^{\top}Q^{\top})_{i,j}&=& \frac{1}{\sqrt{k_h}}q^{\top}_{hj}= \frac{c^{\rm out}_{hj}}{\sqrt{k_j}}. \end{eqnarray} (A.3) We further note that $$(DS^{\top})_{ih}= (S^{\top}\overline{D})_{ih}= \frac{1}{\sqrt{k_h}}\delta_h$$, if $$i \in V_h$$, otherwise $$(DS^{\top})_{ih}=0$$. Thus the statement holds. ii. By using the result in i., the proof in [15, Theorem 2.2] applies. □ A.3 Proof of Proposition 4.3 We first need some technical facts to prove Proposition 4.3. Proposition A.1 Let $$A$$ be an $$n \times n$$ irreducible and non-negative matrix and let $$D=\mathop{\rm diag}(\delta_1, \,{\ldots}\,,\delta_n)$$. Then it holds: i. $$A-D$$ is irreducible, for each $$(\delta_1, \,{\ldots}\,,\delta_n)$$. ii. There exists an eigenvector $$w$$ of $$A-D$$ such that $$w > 0$$ (i.e. each component $$w_i>0$$, $$i=1, \dots, n$$) and the corresponding eigenvalue is $$r(A-D)$$, for each $$(\delta_1, \,{\ldots}\,,\delta_n)$$. Proof. i. From [28]: a $$n \times n$$ matrix $$A$$ is said to be irreducible if for any proper subset $$S \subseteq \left\{1, \ldots, n\right\}$$ there exists $$i \in S$$ and $$j \in S'=\left\{1, \ldots, n\right\}-S$$ such that $$a_{ij} \neq 0$$; since $$A$$ is irreducible, the definition applies immediately to $$A-D$$; ii. See [28, Lemma 4.2]. □ With these background statements we prove Proposition 4.3. Proof. Basically, by Theorem 2.1, we have to show that \begin{equation}\label{eigA} r(A^{\top} - D)=r(\tilde{Q}-\overline D)=r(Q^{\top}- \overline D). \end{equation} (A.4) We first prove that \begin{equation}\label{eqeigA} r( Q^{\top}-\overline D)=r( A^{\top} - D). \end{equation} (A.5) Let $$c \in \mathbb{R}$$ such that both $$a^{\top}_{zz}- \delta_z+ c \geq 0$$, for all $$z=1, \ldots, N$$ and $$q^{\top}_{ii} -\overline \delta_i +c\geq 0$$ for all $$i=1 \ldots, n$$. Let us define $$A^{\top}-D^{\top}+cI_{N \times N}=\hat{A^{\top}}$$ and $$Q^{\top}- \overline{D} + cI_{n \times n}=\hat{ Q^{\top}}$$. $$\hat{A^{\top}}$$ and $$\hat{Q^{\top}}$$ are non-negative and irreducible matrices (see i. in Proposition A.1). We order the eigenvalues of $$\hat{ Q^{\top}}$$ so that $$|\lambda_1(\hat{ Q^{\top}})|\geq |\lambda_2(\hat{ Q^{\top}})| \geq \ldots \geq|\lambda_n(\hat{ Q^{\top}})|$$, and similarly for $$\hat{ A^{\top}}$$. By the Perron–Frobenius theorem [55, Chapter 8], the eigenvalue of maximum modulus of an irreducible and non-negative matrix is real and positive and its corresponding eigenvector, the Perron vector, is the unique (up to a factor) strictly positive eigenvector of the matrix. Hence there exists an eigenvector $$w > 0$$ of $$\hat{ Q^{\top}}$$ corresponding to $$\lambda_1(\hat{ Q^{\top}})$$, that is $$w_i>0$$, for all $$i=1, \,{\ldots}\,, n$$. Now, by ii. in Lemma 4.1 and since $$S^{\top} I_{n \times n}=I_{N \times N} S^{\top}$$, we have that $$S^{\top} w>0$$ is the eigenvector of $$\hat{ A^{\top}}$$ corresponding to $$\lambda_1(\hat{ Q^{\top}})$$. However, because $$S^{\top} w$$ is strictly positive, it must be the Perron vector of $$\hat{ A^{\top}}$$, consequently $$\lambda_1(\hat{A^{\top}})=\lambda_1(\hat{ Q^{\top}})$$. Since $$ r(\hat{ Q^{\top}})=\lambda_1(\hat{ Q^{\top}})= \lambda_1(\hat{ A^{\top}})=r(\hat{A^{\top}}), $$ and \begin{equation}\label{real} r(Q^{\top}-\overline D)+c=r(\hat{ Q^{\top}})=r(\hat{ A^{\top}})=r(A^{\top}-D) + c, \end{equation} (A.6) (A.5) is proved. Now we prove that \begin{equation}\label{eqeigQ} r(\tilde{Q}-\overline D)=r(Q^{\top}-\overline D). \end{equation} (A.7) Let the matrix $$\Lambda=\mathop{\rm diag}(k_i)$$, $$i=1, \ldots, n$$. For any $$n$$-dimensional vector $$v$$ and scalar $$\lambda \in \mathbb{C}$$ we have that \begin{gather*} \left(\tilde{Q} - \overline D\right)v = \lambda v \Leftrightarrow \left( \Lambda^{-\frac12}Q \Lambda^{\frac12}-\overline D\right)v=\lambda v \Leftrightarrow\\ \left(Q \Lambda^{\frac12} - \overline D \Lambda^{\frac12} \right)v = \lambda \Lambda^{\frac12}v\Leftrightarrow\\ \left(Q -\overline D\right) \left(\Lambda^{\frac12} v\right) =\lambda \left(\Lambda^{\frac12}v\right)\!, \end{gather*} hence $$\lambda \in \sigma(\tilde{Q} - \overline D) \Longleftrightarrow \lambda \in \sigma( Q - \overline D)$$, so that (A.7) is verified. In conclusion, from (A.7) and (A.5) it follows (A.4). □ A.4 Proof of Proposition 5.4 Proof. i. Let $$\delta_i=0$$ for some $$i=1,\ldots n$$ and assume that $$\lambda_1(A-D) < 0$$: for the vector $${\mathbf e}_i$$ of the canonical basis it holds $${\mathbf e}_i^{\top} {(A - D)}{\mathbf e}_i={\mathbf e}_i^{\top} A {\mathbf e}_i \geq 0$$ which is a contradiction. ii. The eigenvalues of such kind of matrix vary with continuity with the entries of the matrix [55, Appendix D]. iii. Let us consider $$c>0$$ such that $$-d_{ii}+c\geq 0$$ for all $$i=1,\ldots,n$$. Then, $$A-D+cI$$ is non-negative irreducible and $$\lambda_1(A-D+c\cdot I_n)$$ is actually its Perron–Frobenius eigenvalue [55, Chapter 8]. Now, we can write for any $$\epsilon>0$$ \[ \lambda_1(A-D+\epsilon\, {\mathop{\rm diag}({\mathbf e_i})})=\lambda_1(A-D+\epsilon \, {\mathop{\rm diag}({\mathbf e_i})}+c\, I_n)-c > \lambda_1(A-D +c\,I_n)-c=\lambda_1(A-D), \] where the strict majorization holds because the Perron–Frobenius eigenvalue is strictly increasing in any entry of the matrix [55, 56]. □ A.5 Examples 1. A simple example of the optimal solution for the case of a community network is that of a star graph, where we have two communities, one formed by the central node and the other by the leaf nodes. Assuming that the infection rate is $$\beta$$, we have to find the value of $$\delta_0$$ and $$\delta_1$$ for which $$ \beta Q- D$$ has the maximal eigenvalue which is equal to zero. The characteristic polynomial of $$\beta Q- D$$ for a star graph with $$k$$ leaves is \[ p_{\lambda}(\beta Q- D)=\lambda^2 + (\delta_0 + \delta_1)\lambda +\delta_0\delta_1 - \beta^2 k. \] We observe that $$\lambda=0$$ belongs to the spectrum of $$\beta Q- D$$ if and only if $$\delta_0=\beta^2 k/\delta_1$$. This also ensures that the second eigenvalue is negative and, consequently, $$\lambda=0$$ must be the largest eigenvalue of $$\beta Q- D$$. Thus replacing the value of $$\delta_0$$ obtained above in $$U(\delta_0,\delta_1)= c_0 \delta_0 + c_1 k \delta_1,$$ and setting to zero the following derivative $$U'(\delta_1)=-\frac{c_0 \beta^2 k^3}{\delta_1^2}+ c_1 k,$$ we have that the linear cost optimization is solved for \begin{equation}\label{eq:immtwolev} \delta_0 = \beta k \sqrt{\frac{c_1}{c_0}}, \quad \delta_1= \beta \sqrt{\frac{c_0 }{c_1}}, \end{equation} (A.8) which in turn provides the optimal cost \begin{equation*}\label{eq:starcost} U^*=\beta k \left( c_0\sqrt{\frac{c_1}{c_0}} + c_1\sqrt{\frac{c_0}{c_1}}\right)= 2 \beta k \left(\sqrt{c_1 c_0}\right)\!. \end{equation*} We observe that the optimal cost is linear in the terminal community size $$k$$. In the case of a uniform curing policy, where all nodes are cured at rate $$\delta$$, we have that the value of $$\delta$$ such that $$\lambda=0$$ is the largest eigenvalue of $$\beta Q - D$$ is equal to $$\beta \sqrt{k}$$, thus the value $$U_u$$ of the total cost is \[ U_u = \beta \sqrt{k} (c_0 + c_1 k). \] It is easy to see that the ratio $$U_u/U^*$$ is increasing in $$(1, \infty)$$, moreover we can observe that \begin{equation}\label{eq:ratio} \frac{U_u}{U^*} = O(\sqrt k). \end{equation} (A.9) Hence it is clear that we have an advantage in considering a two-level curing policy, with respect to the uniform curing policy. 2. Now we consider an interconnected star network with two linked central nodes, where each terminal community has the same number of elements $$k$$. We set $$\beta$$ as the infection rate between the central nodes and $$\varepsilon \beta$$ the infection rate between a central node and a node in its adjacent terminal community, where $$\varepsilon > 0$$. After computing the characteristic polynomial of $$\beta Q- D$$, we can see that the zero eigenvalue belongs to the spectrum of $$\beta Q- D$$ provided that $$\delta^2_0 \delta^2_1 - 2 \beta^2 \delta_0 \delta_1 \varepsilon^2 k + \varepsilon^4 \beta^4 k^2- \beta^2 \delta_1^2 =0.$$ The values of $$\delta_0$$ for which $$\lambda=0$$ corresponds to the largest eigenvalue of $$\beta Q- D$$ is equal to $$\delta_0 = \frac{\beta^2 \varepsilon^2 k}{\delta_1} +\beta$$ and the linear cost optimization is solved for \[ \delta_0 = \beta \left(\varepsilon k \sqrt{\frac{c_1}{c_0}}+1\right)\!, \qquad \delta_1=\varepsilon \beta \sqrt{\frac{c_0}{c_1}}. \] Consequently the optimal cost is \begin{equation}\label{twostarcost} U^*= 2 \beta c_0 \left(\varepsilon k \sqrt{\frac{c_1}{c_0}} + 1 \right) + 2 c_1 k \sqrt{\frac{c_0}{c_1}} \varepsilon \beta = 2 \beta \sqrt{c_0c_1}(\varepsilon (k+1) + c_0). \end{equation} (A.10) In the case of a uniform curing policy we have that the value of $$\delta$$ such that $$\lambda=0$$ is the largest eigenvalue is $$(\beta+\sqrt{\beta^2 + 4 \beta^2 \varepsilon^2 k})/2$$ and the value of the total cost is \[ U_u=\hspace{-0.2 em}c_0\hspace{-0.2 em}\left(\hspace{-0.2 em}\beta\hspace{-0.2 em} +\hspace{-0.2 em} \sqrt{\beta^2\hspace{-0.2 em} +\hspace{-0.2 em} 4\beta^2 \varepsilon^2 k}\hspace{-0.2 em}\right)\hspace{-0.2 em} +\hspace{-0.2 em} c_1 k \left(\hspace{-0.2 em} \beta\hspace{-0.2 em} +\hspace{-0.2 em} \sqrt{\beta^2\hspace{-0.2 em} +\hspace{-0.2 em} 4\beta^2 \varepsilon^2 k}\hspace{-0.2 em}\right)\!. \] The ratio $$U_u/U^*$$ is increasing in $$(0,\infty)$$, and again we have that \[ \frac{U_u}{U^*} = O(\sqrt k). \] 3. Now we consider a complete bipartite graph. Basically, we have two communities and we denote by $$k_0$$ the number of elements in the community whose nodes have recovery rate $$\delta_0$$ and $$k_1$$ the number of elements in the community whose nodes has recovery rate $$\delta_1$$. The optimal curing rates are \[ \delta_0= \beta k_1 \sqrt{\frac{c_1}{c_0}}, \qquad \delta_1= \beta k_0 \sqrt{\frac{c_0}{c_1}}, \] and the optimal cost is \[ U^*= c_0 k_0 \beta k_1 \sqrt{\frac{c_1}{c_0}} + c_1 k_1 \beta k_0 \sqrt{\frac{c_0}{c_1}}. \] In the case of a uniform curing policy the value of $$\delta$$ such that $$\lambda=0$$ is the largest eigenvalue is \[ \delta= \beta \sqrt{k_0 k_1}, \] and the cost is \[ U_u= c_0 k_0 \beta \sqrt{k_0 k_1}+ c_1 k_1 \beta \sqrt{k_0 k_1}. \] In this case the asymptotic behaviour of $$U_U/U^*$$ for high values of $$k_0$$ and $$k_1$$ depends on the direction in which we move, thus we cannot say anything in this regard. © The authors 2017. Published by Oxford University Press. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Optimal curing policy for epidemic spreading over a community network with heterogeneous population
JF - Journal of Complex Networks
DO - 10.1093/comnet/cnx060
DA - 2018-10-01
UR - https://www.deepdyve.com/lp/oxford-university-press/optimal-curing-policy-for-epidemic-spreading-over-a-community-network-0LuX0xJYWy
SP - 800
VL - 6
IS - 5
DP - DeepDyve
ER -