A GEOMETRIC PROGRAMMING APPROACH FOR OPTIMAL RESOURCE ALLOCATION TO CONTROL EPIDEMIC OUTBREAKS IN ARBITRARY NETWORKS

This paper proposed a control-strategies for nodes to control the spread of an epidemic outbreak in arbitrary directed graphs by optimally allocating their resources throughout the network. Epidemic propagation is well modeled as a networked version of the Susceptible-Exposed-Infected-Susceptible (SEIS) epidemic process. Using the Kolmogorov forward equations and mean-field approximation, we present a mean-field model to describe the spreading dynamics and prove the existence of a necessary and sufficient condition for global exponential stability. Based on this stability condition, we can derive another condition to control the spread of an epidemic outbreak in terms of the eigenvalues of a matrix that depends on the network structure and the parameters of the model. According to different control purposes and conditions, two types of control-theoretic decision can be considered: 1)given a fixed budget, find the optimal resource allocation to achieve the highest level of containment, 2)given a decay rate of epidemic, find the minimum cost to control the spreading process at a desired decay rate. A geometric program can be formulated to solve the optimal problems and the existence of solutions is also proved. Numerical simulations can illustrated our results.


Introduction
Development of strategies to control spreading processes in networks has brought much attention due to its applications in many relevant fields, including computer virus [1], public health [2][3], and information propagation over social network [4].The dynamic of spreading processes in networks not only depend on the epidemic model but also the structure of the contact network.In this context, the spread process is modeled by a variant of SIS epidemic model that includes a state of "expose".The individual infected rate and recovery rate can be modified within a feasible range by allocating resources in each node.Based on this model, an efficient convex framework can be formulated to solve the optimal problems.The dynamic behavior of spreading processes in networks has been widely studied.Meanfield models over arbitrary contact graphs were brought to the forefront in both continuous time [5] and discrete-time.[6] is a paper about the dynamic behavior of spreading processes in arbitrary contact networks for the case of discrete-time dynamic.[7] consider a continuous time SIS model over arbitrary graphs using mean field theory and provide a condition for globally asymptotically stable of the disease-free state.It should be noted that most models are developed base on undirected graphs, however, in practice, directed graphs may be more appropriate to the spread of diseases in human populations.Therefore, we study the spread processes in an arbitrary directed network of heterogeneous nodes.
Designing control strategies of spread processes in networks is a significant problem.Several papers have proposed some methods on different aspects of this problem.In [8], the authors proposed a semidefinite programming (SDP) to find the optimal strategies of resource allocation in an undirected network.[9] uses a linear-fractional optimization program to find the optimal investment on disease awareness in a social network.In [10], Borgs et al. provides a probabilistic analysis for the case of a given contact network to characterize the optimal allocation of a fixed amount of antidote.[11] provided a eigenvalue sensitivity analysis ideas to design optimal strategies of allocation resource to control the spread of a virus.The relationship between the recovery parameters and distributed approach is explored in [12].Our work is based on [13] and [14], in which a continuous-time time Markov processes, called the N-intertwined model, is used to analyze and control the spread of a SEIV epidemic model.This paper's layout is as follows.In section 2, some important notation and background need to be introduced.In section 3, we formulate two resource allocation problems for epidemic propagation in the network.In section 4, a convex optimization framework is used to efficiently solve the allocation problems both on strongly connected graphs and general directed graphs (not necessarily strongly connected).In section 5, simulations illustrate our results.

Notation and preliminaries
Some graph-theoretical nomenclature and the dynamic spreading model need to be introduced in this section.

Graph Theory
Let G = (V , E ) denote an directed graph, where V = {v 1 , •, v n } denote the set of n nodes and E ⊆ V × V denote the set of ordered pairs of nodes called directed edges.An edge from e j pointing towards e i denote as (e j , e i ).Let N in i = { j : (e j , e i ) ∈ E } denote the in-neighborhood of node i.A directed graph is strongly connected if and only if there is a directed path from e j to e i for every pair of nodes e i , e j ∈ V .
The adjacency matrix of a digraph G , denoted by U G = [u i j ], is an n × n matrix defined entry-wise as u i j = 1 if edge (e i , e j ) ∈ E , and u i j = 0 otherwise.Hence, the adjacency matrix U G = [u i j ] is always nonnegative.An nonnegative adjacency matrix is irreducible if and only if its associated graph is strongly connected.Given an n × n matric Z, the sets of eigenvalues and corresponding eigenvectors of Z are defined as λ 1 (Z), •, λ n (Z) and η 1 (Z), •, η n (Z), where we order eigenvalues in decreasing order of their real parts, i.e., R(λ λ 1 (Z) and η 1 (Z) are called the dominant eigenvalue and eigenvector of Z. ρ(Z) is defined as the spectral radius of Z and equal to the maximum modulus across all eigenvalues of Z.
It should be note that we only consider unweighted digraphs in this paper, hence, the adjacency matrix is always nonnegative.

Stochastic heterogeneous SEIS model
The Susceptible-Exposed-Infected-Susceptible epidemic model is a variant of the SIS model by including a state of "exposed" which the node has been exposed to the disease and is contagious, but is not aware of the contagion.This model is a continuous-time markov process and each node in the network can be in one out of three possible states.Consider a network of n individuals described by the adjacency matrix U G = [u i j ], and the parameters are defined as follows: β i : infection rate that the susceptible node i transitions to the exposed state by contact with its exposed neighbors.α i : infection rate that the susceptible node i become exposed by contact with its infected neighbors.
γ i : the rate at which the exposed node i be infected.
δ i : the recovery rate of node i.
Let v k (τ) ∈ {0, 1} and v j (τ) ∈ {0, 1}.Where v k (τ) = 1 indicates that node k is in the infected sate and other states of node k at time τ denote as v k (τ) = 0. Similarly, v j (τ) = 1 indicates that node j is in the exposed sate, and v j (τ) = 0 otherwise.Three possible types of stochastic transitions during the time interval [τ + ∆τ): a): Assuming node i is in the susceptible state at time τ.This node can switch to the exposed state during the small time interval [τ, τ + ∆τ) with a probability that depends on: (i) its infection rates β i and α i ; (ii) the strength of its incoming connections {u i j , for j ∈ N in i }; (iii) the states of its in-neighbors {v j (τ), for j ∈ N in i } and {v k (τ), for k ∈ N in i }.Formally, the probability of this transition is given by where ∆τ > 0 is an arbitrarily small time interval, and b): Assuming node i in exposed state, the probability of i transit to infected in the time interval [τ, τ + ∆τ) is given by c): Assuming node i is infected, the probability of i recovering back to the susceptible state in the time interval [τ, τ + ∆τ) is given by The Markov process with 3 n states described in the above is very hard to analyze due to the exponential size of the state space.Therefore, we use a mean-field approximation of its dynamics.This approximation is widely used in the field of epidemic analysis and control, since it performs numerically well for many realistic network topologies.
Let h i and g i denote the probabilities of node i to be exposed and infected, respectively.
Then, the probabilities of node i to be susceptible is 1 − h i − g i .Using the Kolmogorov forward equations and a mean-field approach, one can approximate the dynamics of the epidemic spread using a system of 2n ordinary differential equations, as follows: ) We can write the mean-field approximation equations in matrix form as where Proposition 1 Consider the heterogeneous SEIS epidemic model in (2.3) and assume that U G ≥ 0 and for some ς > 0, the disease-free equilibrium is globally exponentially stable, and ς is called exponential decay rate.

Problem formulation
To control the spread of an epidemic in a given network, the work of proposing an efficient optimization framework to find the optimal resource distribution is very significant.In this paper, we consider two types of resources: 1) preventive resource (e.g.vaccines to reduce infection rate β i , α i ), 2) corrective resources (e.g.antidotes to help increase recovery rates δ i ).
We can distribute resources to modify the parameters β , α, δ , and γ, within the feasible ranges as follows: where β i and α i are the minimum possible infection rates for node i.They can be achieved by allocating a large amount of vaccines at node i. βi and ᾱi are the maximal infection rates without any preventive resource for node i.Similarly, the δ i is a natural recovery rate of node i, and δi is the maximal recovery rate which is achieved with enough corrective resources.For convenience, let ∇ = {β , α, γ, δ } denote as the global set of parameters.

The cost of preventive and corrective resources
We define three cost functions, the vaccination cost function f i (β i ), g i (α i ) and treatment cost function h i (δ i ).The cost functions are node-dependent and present the following properties: (1) Assuming the vaccination cost f i (β i ) and g i (α i ) are monotonically decreasing in the interval [β i , βi ], [α i , ᾱi ], respectively.Antidote cost g i (δ i ) is monotonically increasing with regard to δ i .
(2) In the absence of investment, Apart from the above properties, we assume the cost functions have the following forms to obtain a tractable convex framework.
Note that we have normalized these cost functions to have values in the interval [0, 1].
Therefore, the cost functions of preventive resource, f i (β i ), g i (α i ) and the corrective resource cost, h i (δ i ), are twice differentiable and satisfies the following constrain: Notice that, since f i , g i are monotonically decreasing, h i is monotonically increasing, we have that f i < 0, g i < 0 and h i > 0. The results implies that f i > 0, g i > 0 and h i < 0. Therefore, the assumption is stronger than convexity.

Problem statements
In this section we present two types of resource allocation problems for SEIS model: (1) the budget-constrained allocation problem.Our aim is finding the optimal allocation of vaccination and antidotes to maximize the exponential decay rate ς when given the total budget T > 0.
(2) the rate-constrained allocation problem.Given the exponential decay rate ς , find the cost-optimal distribution of vaccines and antidotes to eradicate the disease with decay rate greater than equal to ς .
Problem 1(Budget-constrained allocation) Given the total budget T , this problem can be state as the following optimization problem: for all i = 1, •, n, where (3.3) is the budget constraint.
Problem 2 (Rate-constrained allocation) Given a desire decay rate ς , the rate-constrained allocation is formulated as follows: for all i = 1, •, n, where (3.8) constrains the decay rate to ς .
In the following section we propose an approach to solve these problems in polynomial time.

4.A convex framework for optimal resource allocation
A convex formulation can be use to solve both the rate-constrained allocation problem and the budget-constrained problem in unweighted, directed networks using GP.We first solve problems in the strongly connected digraphs, then extend the results to general digraphs.
Some important concepts need to be briefly reviewed.Denote n decision variables ξ = (ξ 1 , •, ξ n ), ξ i > 0. In the context of geometric programs, a monomial function h(ξ ) is defined as a real-valued function of the form h(ξ , where a > 0 and k i ∈ R for all i = 1, •, n.The sum of monomials is defined as the polynomial function, i.e., q(ξ n with a i > 0 and k j,i ∈ R for all j ∈ {1, •, n} and i ∈ {1, •, I}.
A geometric problem is an optimization problem of the form: where f (ξ ) and q i (ξ ) are polynomial functions, h i (ξ ) are monomials.
Note that polynomials and monomials are convex in log-scale, therefore f is a convex function in log-scale.The quasiconvex optimization problem GP can be transformed to a convex problem based on a logarithmic change of variable ϕ i = logξ i , and a logarithmic transformation of the objective and constraint functions.After this transformation, the GP in (4.1) takes the form as follows: minimize F(ϕ) where F(ϕ) = log f (exp ϕ) and Q i (ϕ) = log q i (exp ϕ).Also, assuming that n , the equality constraint after the logarithmic change in variable can be obtained In summary, (4.2) is a convex optimization problem in standard form and can be efficiently solved in polynomial time.
In the following section, we show how to transform our problems into GPs.In our transformation, the theory of nonnegative matrices and the Perron-Frobenius lemma are very useful.

GP for strongly connected digraphs
For a strongly connected digraph J , its adjacency matrix Z is irreducible.Then the follow lemma holds for the spectral radius of the adjacency matrix of any unweighted, strongly connected digraph.
Lemma 1 (Perron-Frobenius):Z is a nonnegative, irreducible matrix.Then, the following statements about its spectral radius ρ(Z) hold: (1) ρ(Z) > 0 is a simple eigenvalue of Z, (2) Zω = ρ(Z)ω, for some ω > 0, Based on the above results, we have the following result: Proposition 2: Consider the n × n nonnegative, irreducible matrix Z(x) with entries being either polynomials or 0 with domain ξ ∈ Ω, where Ω is defined as polynomials.Then, we can solve the following GP via minimize λ 1 (Z(ξ )): Based on the above results, assuming that the contact graph J is strongly connected, we can solve both the budget-constrained and the rate-constrained problems.
Theorem 1(Solution to the budget-constrained problem): For strongly connected digraphs, problem 2 can be solved by the following GP min ) First, note that the matrix L is not a nonnegative, we can define a nonnegative matrix from L by simply adding a constant k = max{γ i , δi Notice that the matrix L is nonnegative and k = max{γ i , δi } for i = 1, •, n.Then, according to proposition 2, we have that maximizing ς in (3.1) is equivalent to minimizing λ 1 ( L) under budget constraints.We can minimizing λ 1 ( L) by solving the following GP: for all i ∈ {1, •, n}.The matrix L is nonnegative and irreducible if the adjacency matrix U J corresponding to a strongly connected digraph.Therefore, applying proposition 2, the first constraint can be rewrite as the following two constraints: for all i ∈ {1, •, n}, where δi = k − δ i .After the above change of variables, the problem (2) rewrite as a standard GP form (4.1).We can find the optimal resources allocation in a heterogeneous network under the budget constraint.
Theorem 2(Solution to the rate-constrained problem): problem 2 can be solved by solving the following GP: for all i ∈ {1, •, n}, where δ Proof: The proof is similar to the one for theorem 1, so we omit it.
In this section, we have provided the solutions to problems for strongly digraphs.Then we show how to extend it to general connected digraphs.

GP for general connected digraphs
The adjacency matrix U is irreducible if and only if its graph is strongly contact.So the Perron-Frobenius lemma is not applicable to digraphs that are not strongly connected.For general digraphs, the statements in the P-F lemma are weaken, as follows: Lemma 2 (Perron-Frobenius):Z is a n × n nonnegative matrix.Then, the following statements about its spectral radius ρ(Z) hold: (1) ρ(Z) ≥ 0 is an eigenvalue of Z, (2) Zω = ρ(Z)ω, for some ω ≥ 0, Note that the components of ω in proposition 2 are strictly positive, however, the components of ω in lemma 2 are nonnegative.So if we want to use GP, this issue must be resolved.Defined the sets Θ = {i : η i = 0 and ω i = 0}.Hence the variables for i / ∈ Θ can be excluded from the GP's in theorems 1 and 2. Hence the allocation problems can be split into two different sets of decision variables.
Rate-Constrained Allocation Problem for General Digraphs: If i / ∈ Θ, the optimization problem holds for the variables: Thus, for f i , g i decreasing and h i increasing, it is obviously that the minimum investment for all i / ∈ Θ correspond to the optimal infection rates βi , ᾱi and optimal recovery rate δ i .
Theorem 3: On the other hand, for i ∈ Θ, the optimal solution can be obtained from the following GP: ) Theorem 4: Budget-Constrained Allocation Problem for General Digraphs: For i / ∈ Θ, it's easily to find the optimal spreading and recovery rates are βi , ᾱi and δ i .Since one of the nodes with zero eigenvalue.Thus, we can rewrite the budget-constrained allocation problem as following GP for general digraphs: for i ∈ Θ, where δ * i = k − δi * and λ * 1 (L) ≤ λ * − k.Theorem 3 and 4 provided the solutions to both two problems for general digraphs.

Simulation
We have developed an optimization program for determining optimal-cost parameter distributions such that the desired equilibrium is stabilized.In this section, we present how the geometric programming solve the optimization problem of resource allocation by simulating a In Fig. 1 we plot the cost functions which are given in subsection 3.1.Here the abscissa is the amount investment on the node n i and the ordinates are the infection (red and blue line) and recovery (black line) rates achieved by the investment.As we increase the amount invested on protection resources from 0 to 1, the infection rate of that node is reduced from ( βi , ᾱi ) to (β i , α i ) (red and line).Similarly, as we increase the amount invested on corrective resources at a node n i , the recovery rate grows from δ i to δi (black line).Fig. 2 demonstrates the performance achieved by using the rate-constrained allocation algorithm from Theorem 2. We have that in the optimal allocation some nodes receive no resources at all; some nodes receive only preventive or corrected resources, and some nodes receive a mixture of preventive and corrective resources.

Conclusion
We have presented a convex optimization framework to find the optimal allocation resources of the SEIS epidemic model on arbitrary directed graphs.A necessary and sufficient condition for global exponential stability can be derived from the eigenvalues of a matrix that depends on parameters of the model and the network structure.Furthermore, We have formulated optimization programs for determining optimal resources allocation and reformulating them as geometric programs that can efficiently solve the optimal resource allocation problem.For future work we plan to study the endemic equilibrium which the disease-free equilibrium is not globally asymptotically stable.

FIGURE 1 .
FIGURE 1.A plot of infection rate β i (in red), α i (in blue) and recovery rate δ i (in black) achieved at node n i after the investments on preventive and corrective resources are made on that node.

FIGURE 2 .
FIGURE 2. The optimal investment on preventive and corrective resources for ten nodes in a strongly connected digraph, where the abscissas are nodes and the ordinates are the amount invested on preventive and corrective resources.