Optimal networks for dynamical spreading
aa r X i v : . [ phy s i c s . s o c - ph ] J a n Optimal networks for dynamical spreading
Liming Pan, Wei Wang, ∗ Lixin Tian, and Ying-Cheng Lai
4, 5 School of Computer and Electronic Information, Nanjing Normal University, Nanjing, Jiangsu, 210023, China Cybersecurity Research Institute, Sichuan University, Chengdu 610065, China School of Mathematical Sciences, Nanjing Normal University, Nanjing, Jiangsu 210023, China School of Electrical, Computer and Energy Engineering,Arizona State University, Tempe, Arizona 85287, USA Department of Physics, Arizona State University, Tempe, Arizona 85287, USA (Dated: January 7, 2021)The inverse problem of finding the optimal network structure for a specific type of dynamical processstands out as one of the most challenging problems in network science. Focusing on the susceptible-infected-susceptible type of dynamics on annealed networks whose structures are fully characterized by the degree distri-bution, we develop an analytic framework to solve the inverse problem. We find that, for relatively low or highinfection rates, the optimal degree distribution is unique, which consists of no more than two distinct nodal de-grees. For intermediate infection rates, the optimal degree distribution is multitudinous and can have a broadersupport. We also find that, in general, the heterogeneity of the optimal networks decreases with the infectionrate. A surprising phenomenon is the existence of a specific value of the infection rate for which any degreedistribution would be optimal in generating maximum spreading prevalence. The analytic framework and thefindings provide insights into the interplay between network structure and dynamical processes with practicalimplications.
I. INTRODUCTION
In the study of dynamics on complex networks, most previ-ous efforts were focused on the forward problem: How doesthe network structure affect the dynamical processes on thenetwork? The approaches undertaken to address this ques-tion have been standard and relatively straightforward: Oneimplements the dynamical process of interest on a given net-work structure and then studies how alterations in the networkstructure affect the dynamics. The dynamical inverse prob-lem is much harder: finding a global network structure thatoptimizes a given type of dynamical processes. Despite theextensive and intensive efforts in the past that have resultedin an essential understanding of the interplay between dynam-ical processes and network structure, previous studies of theinverse problem were sporadic and limited to a perturbationtype of analysis, generating solutions that are at most locallyoptimal only [1, 2]. The purpose of this paper is to present anddemonstrate an analytic framework to address the dynamicalinverse problem.To be concrete, we will focus on spreading dynamics onnetworks for which a large body of literature has been gen-erated in the past on the forward problem, i.e., how networktopology affects the characteristics of the spreading, such asthe outbreak threshold and prevalence [3, 4]. For example,under the annealed assumption that all nodes with the samedegree are statistically equivalent, it was found [5] that the epi-demic threshold of the susceptible-infected-susceptible (SIS)process is given by h k i / h k i , where h k i and h k i are the firstand second moments of the degree distribution, respectively.In situations where the second moment diverges, the thresh-old value is essentially zero, meaning that the presence of a ∗ [email protected] few hub nodes can greatly facilitate the occurrence of an epi-demic outbreak. An understanding of the interplay betweenthe network structure and the spreading dynamics is essentialto articulating control strategies. For example, the importantrole played by the hub nodes suggests a mitigation strategy:Vaccinating these nodes can block or even stop the spread ofthe disease [6, 7]. Likewise, if the goal is to promote infor-mation spreading, then choosing the hub nodes as the initialseeds can be effective [8, 9].The inverse problem is motivated by the application scenar-ios in which one strives to optimize the network structure toachieve desired or improved performance [10]. Optimizationand invention have been applied to problems such as virusmarketing [11], social robots detection [12], containment offalse news spreading [13], and polarization reduction in so-cial networks [14]. For spreading dynamics on networks, thefew existing studies are focused on applying small perturba-tions to the network structure to modulate the dynamical pro-cess [1, 2]. From the point of view of optimization, since theperturbations are local, the resulting solution is locally opti-mal at best.We address the following questions: Does a globally op-timal network exist and if yes, can it be found to maximizethe prevalence of the spreading dynamics? Such a networkis necessarily extremum. For general types of spreading dy-namics, to analytically solve this inverse problem is currentlynot feasible. However, we find that the SIS type of spreadingdynamics does permit an analytic solution. In particular, theannealed approximation stipulates that the network structurecan be fully captured or characterized by its degree distribu-tion. The problem of finding the optimal networks can thenbe formulated as one to find the optimal degree distributionthat maximizes the prevalence of the SIS spreading dynamics,which can be analytically solved by exploiting the heteroge-neous mean-field (HMF) theory [3]. Notwithstanding the ne-cessity of imposing the annealed approximation to enable an-alytic solutions, the essential physical ingredients of the SISdynamics are retained.Our main results are the following. Taking a variational ap-proach to solving the HMF equation, we obtain a necessarycondition for the optimal degree distribution. The conditiondefines a set of candidate optimal degree distributions, and weshow that a degree distribution is globally optimal if and onlyif it belongs to the set. However, if the set is empty, which canoccur for relatively low and high infection rates, the necessarycondition stipulates that a local extremum distribution mustconcentrate on no more than two distinct nodal degree valuesthereby substantially narrowing the search for the optimal net-work. Searching through all possible distributions under theconstraint leads to the optimal degree distribution that can beproved to be unique. For intermediate infection rates, multi-ple optimal degree distributions with a broader support exist,which lead to identical spreading prevalence. In addition, ourtheory predicts the existence of a particular value of the in-fection rate for which every degree distribution is optimal. Ageneral trend is that the degree heterogeneity of the optimaldistribution decreases with the infection rate.Our paper represents a first step toward finding a global op-timal network structure for spreading dynamics. From a the-oretical point of view, developing a method to find such ex-tremum networks represents a feat that would provide deeperinsights into the interplay between network topology andspreading dynamics. From a practical perspective, the solu-tion can be exploited to design networks that are capable ofspreading information or transporting material substances inthe most efficient way possible.In Sec. II, we introduce the HMF theory for the SIS dy-namics and set up the basic framework for the optimizationproblem. In Sec. III, we employ a variational method to de-rive the necessary condition for a degree distribution to be anextremum among all feasible distributions. Solutions of theoptimal degree distribution are presented in Sec. IV, and itsproperties are discussed in Sec. V. The paper is concluded inSec. VI with a discussion. II. PROBLEM FORMULATION AND SKETCH OF MAJORMATHEMATICAL STEPS
In the SIS model, each node can be either in the suscep-tible or in the infected state, and we assume the nodal stateevolves continuously with time. During the spreading pro-cess, a susceptible node is infected by its neighbors with therate λ , whereas an infected node recovers at the rate γ . Tostudy the equilibrium properties of the dynamical process, itis convenient to set γ = 1 so that λ is the sole dynamicalparameter.In the HMF theory, all the nodes with the same degree arestatistically equivalent [3]. Consider a vector of nodal degrees k ≡ [ k , k , · · · , k n ] T , where the elements are arranged in adescending order: k > k > · · · > k n . The degree dis-tribution is fully specified by a probability vector defined as p ≡ [ p , p , · · · , p n ] T , where p i ≥ is the probability that arandomly chosen node has degree k i . Let x i ( t ) be the prob- ability that a node with degree k i is infected at time t . Giventhe probability vector p , the HMF equation is dx i ( t ) dt = − x i ( t ) + λk i [1 − x i ( t )] Θ (1)for i ∈ { , · · · , n } , where Θ = 1 h k i n X j =1 p j k j x j ( t ) . (2)In Ref. [15], it was proved that the HMF equation has aunique global stable equilibrium point x ∗ . In addition, for λ < h k i / h k i , we have x ∗ i = 0 for all i ∈ { , · · · , n } , whereas for λ > h k i / h k i , we have < x ∗ i < for all i ∈ { , · · · , n } .The spreading prevalence in the equilibrium state is ψ ( p ) = n X i =1 p i x ∗ i , (3)where, to simplify the notations, we have omitted the depen-dence of ψ ( p ) on λ . Let P be the family of all degree dis-tributions with a fixed average degree defined on k . That is,with a prespecified constant z > , for any p ∈ P , we have P ni =1 p i k i = z . Our goal is to find p o ∈ P that maximizes ψ ( p ) : p o = argmin p ∈P ψ ( p ) . (4)The optimization problem is nontrivial only when the valueof λ is larger than the epidemic threshold at least for one p ∈ P . The Bhatia-Davis inequality stipulates that the sec-ond moment of p is maximized when p concentrates on theend points k and k n . In this case, the second moment is h k i = zk + zk n − k k n . The optimization problem is non-trivial only when the following condition is met: λ > λ ≡ zzk + zk n − k k n . (5)In this case, if there is a unique solution p such that λ >z/ h k i , it gives the optimal degree distribution p o .Our goal is to analytically find the solutions for the opti-mization problem defined in Eq. (4). As the mathematicalderivations involved are lengthy, it may be useful to sketchthe basic idea, tools used, and the results, which we organizeas the following three major steps.1. Mathematically, Eq. (4) defines a variational problemfor the HMF equations in Eq. (1), which can be studiedthrough the standard calculus-of-variation techniques.In Sec. III A, we adopt a variational approach for theHMF equations in Eq. (1) and derive the necessary con-dition for a degree distribution to be optimal. In par-ticular, we impose a perturbation to the degree distribu-tion as p ′ = p + α ¯ p and derive a formula that predicts ¯ ψ ( p , p ′ ) , the part of the incremental spreading preva-lence which is linear in α . For p to be a candidate max-imum, ¯ ψ ( p , p ′ ) must be nonpositive for any choice of ¯ p , and this leads to the necessary condition for the localminima.2. The next task is to study the necessary condition result-ing from the variational analysis. In Sec. III B, througha sequence of algebraic arguments, we show that forany p satisfying the necessary condition, it is only pos-sible to have either (i) ¯ ψ ( p , p ′ ) = 0 or (ii) ¯ ψ ( p , p ′ ) < for all feasible perturbations. This means that it is im-possible to find a certain p such that ¯ ψ ( p , p ′ ) = 0 and ¯ ψ ( p , p ′′ ) < for p ′ = p ′′ . Further, in Sec. III B, weshow that the condition ¯ ψ ( p , p ′ ) = 0 can be reduced toa linear equation in p [the first equation in (25)] which,together with the probability constraint P ni =1 p i = 1 and the average degree constraint P ni =1 p i k i = z , de-fines a set of candidate optimal degree distributions P o . In Sec. III C, by analyzing the three linear equa-tions, we show that if P o is nonempty, then any p isa global maximum if and only if p ∈ P o . Concur-rently, if P o is empty, the optimal degree distributionwith ¯ ψ ( p , p ′ ) < will concentrate on no more thantwo distinct nodal degrees.3. Finally, in Sec. IV A, we derive the condition when theset P o is nonempty by analyzing the three linear equa-tions defining the set [Eq. (25)]. In particular, P o isnonempty for λ ∈ [ λ , λ ] (see Sec. IV A for explicitdefinitions of λ and λ ). For λ < λ or λ > λ and P o indeed empty, we find the optimal degree distributionsby solving the HMF equations explicitly (Sec. IV B). III. NECESSARY CONDITION FOR LOCAL EXTREMAAND CONSEQUENCES
In this section, we first study the optimization problem de-fined in Eq. (4) using several techniques from the calculus ofvariation. The calculation provides a necessary condition forfinding the local maxima. We then analyze the necessary con-dition in detail to find the global optimal degree distributions.
A. Variational method
We study the variation problem in Eq. (4) using the stan-dard techniques from the calculus of variations. Briefly, weapply a perturbation to the degree distribution p in Eq. (1) andcalculate the linear response for the spreading prevalence. Alocal maximum necessarily has non-positive linear responsesfor any feasible perturbation.For a fixed λ > h k i / h k i , let x ∗ be the corresponding glob-ally stable equilibrium point of the HMF equation. We imposea small variation on p i , p ′ i = p i + α ¯ p i , (6)where ¯ p specifies the direction of the variation and α > con-trols its magnitude. For the perturbed degree distribution to befeasible, i.e., p ′ ∈ P , the following conditions are necessary: n X i =1 ¯ p i = 0 and n X i =1 ¯ p i k i = 0 . (7) In addition, the perturbed degree distribution p ′ must satisfythe probability constraints ≤ p ′ i ≤ .Let x ′ ( t, α ) be the trajectory of the perturbed system. Thetime evolution of x ′ ( t, α ) is described by the HMF equationwith p replaced by p ′ and x i ( t ) in Eq. (1) by x ′ i ( t, α ) . Asshown in Appendix A, x ∗ is a continuously differentiablefunction of p for λ > z/ h k i , enabling the following expan-sion of x ′ ( t, α ) about x ∗ i : x ′ i ( t, α ) = x ∗ i + α ¯ x i ( t ) + o ( α ) , (8)where ¯ x i ( t ) is the response to the perturbation which is linearin α . Taking the derivative with respect to α at α = 0 , we ob-tain ∂x ′ i ( t, α ) /∂α | α =0 = ¯ x i ( t ) . The time derivative of ¯ x i ( t ) is then given by d ¯ x i ( t ) dt = ∂∂α (cid:12)(cid:12)(cid:12)(cid:12) α =0 dx i ( t, α ) dt , (9)which, after some algebraic manipulations, can be rewrittenas d ¯ x ( t ) dt = J ¯ x ( t ) + ξ, (10)where J is the n × n Jacobian matrix that does not dependon ¯ p and ξ is a vector of length n that depends on ¯ p . Theelements of J and ξ are given by J ij = − δ i,j (1 + λk i Θ ∗ ) + λz k i (1 − x ∗ i ) k j p j , (11)and ξ i = λz k i (1 − x ∗ i ) n X j =1 k j ¯ p j x ∗ j , (12)respectively. In Eq. (11), δ i,j is the Kronecker δ and Θ ∗ isobtained by substituting x ( t ) = x ∗ into Eq. (2).Equation (10) defines a linear system with the solution, ¯ x ( t ) = e J t ¯ x (0) + (cid:0) e J t − I (cid:1) J − ξ, (13)where e J t is the matrix exponential of J t . In Appendix B,we show that the eigenvalues of J have negative real parts. Inthe long time limit, we then have ¯ x ∗ = lim t →∞ ¯ x ( t ) = −J − ξ. (14)With the perturbed degree distribution and Eq. (8), we canexpress the spreading prevalence as ψ ( p ′ ) = ψ ( p ) + α ¯ ψ ( p , p ′ ) + o ( α ) , (15)where ¯ ψ ( p , p ′ ) is the part of incremental spreading prevalencethat is linear in α , ¯ ψ ( p , p ′ ) = n X i =1 (¯ p i x ∗ i + p i ¯ x ∗ i ) . (16)Substituting Eq. (14) into Eq. (16), we have (after some alge-braic manipulations) ¯ ψ ( p , p ′ ) = n X i =1 χ i ¯ p i , (17)where χ i is given by χ i = x ∗ i λk i Θ ∗ P nj =1 p j k j (1 − x ∗ j ) P nj =1 p j k j ( x ∗ j ) ! , (18)The detailed derivation of Eq. (18) is presented in Ap-pendix C. The necessary condition for the degree distribu-tion p to be a local maximum is if and only if the inequality ¯ ψ ( p , p ′ ) ≤ holds for all feasible perturbations. B. Consequences of the necessary condition
Equations (17) and (18) allow us to significantly narrow thesearch range for the optimal degree distribution through theprocess of elimination. In the following, we analyze the nec-essary condition by proving that it is only possible to haveeither (i) ¯ ψ ( p , p ′ ) = 0 or (ii) ¯ ψ ( p , p ′ ) < for all feasibleperturbations. That is, it is impossible to find p such that ¯ ψ ( p , p ′ ) = 0 and ¯ ψ ( p , p ′′ ) < for p ′ = p ′′ . We thenshow that ¯ ψ ( p , p ′ ) = 0 can be reduced to an equation thatis linear in p , based on which the spreading prevalence forany p satisfying ¯ ψ ( p , p ′ ) = 0 can be directly obtained with-out solving the HMF equations. The results in this sectionare obtained through algebraic manipulations of the equation ¯ ψ ( p , p ′ ) = 0 .The starting point of our analysis is to determine when thelinear variation ¯ ψ ( p , p ′ ) vanishes. A feasible perturbation ¯ p must satisfy the constraints in (7), so ¯ p must have at least threenonzero elements. Pick any m ≥ points { k i , k i , · · · k i m } from k and consider a perturbation ¯ p whose elements arenonzero only on these points. The linear variation ¯ ψ ( p , p ′ ) vanishes only if Z ( m ) ¯ p = 0 , where Z ( m ) is a × m matrix, Z ( m ) = · · · k i k i · · · k i m χ i χ i · · · χ i m . (19)The first two rows in Z ( m ) correspond to the constraints for ¯ p in (7), while the last row is the result of the definition of ¯ ψ ( p , p ′ ) in Eq. (17). To gain insights, we temporally dis-regard the probability constraint p ′ ∈ [0 1] n (which will beincluded in the analysis later). Under this condition, any ¯ p that makes the linear variation ¯ ψ ( p , p ′ ) vanish belongs to thenull space of Z ( m ) . By the rank-nullity theorem, we have nullity( Z ( m ) ) = m − rank( Z ( m ) ) . The dimension of thespace for all feasible perturbations, i.e., the nullity of the sub-matrix consisting of the first two rows of Z ( m ) , is m − . Asa result, the linear variation vanishes for all directions of per-turbation if nullity( Z ( m ) ) = m − , which further implies thecondition rank( Z ( m ) ) = 2 . We thus have that the linear vari-ation vanishes if and only if the third row of Z ( m ) is a linearcombination of the first two rows. Setting the right-hand-side of Eq. (1) to zero, we obtain theequilibrium solution as x ∗ i = λk i Θ ∗ / (1 + λk i Θ ∗ ) . From thedefinition of χ i in Eq. (18), we have χ i = λk i λk i Θ ∗ λ Θ ∗ k i P nj =1 p j k j (1 − x ∗ j ) P nj =1 p j k j ( x ∗ j ) ! . (20)If the following holds P nj =1 p j k j (1 − x ∗ j ) P nj =1 p j k j ( x ∗ j ) = 1 , (21)then we have χ i = λk i . In this case, the third row of Z ( m ) is exactly the second row multiplying by λ and wehave rank( Z ( m ) ) = 2 . Moreover, if Eq. (21) holds, then rank( Z ( m ) ) = 2 holds for any choice of perturbation with m ≥ . In other words, the linear variation thus vanishes forall directions of perturbation.In the above analysis, we have not required p + α ¯ p ∈ [0 1] n .A direction of perturbation ¯ p would be infeasible if an elementof p has p i = 0 or p i = 1 . Nevertheless, as Eq. (21) guar-antees Z ( m ) ¯ p = 0 for any m , the linear variation ¯ ψ ( p , p ′ ) vanishes in any direction of perturbation, feasible or infea-sible. In fact, in the proof of x ∗ being a continuously dif-ferentiable function of p (Appendix A), it is not necessaryto require p i = 0 or p i = 1 for any i ∈ { , · · · , n } . Thismeans that the perturbation in an infeasible direction can stillbe well-defined, although it is physically irrelevant. Conse-quently, Eq. (21) provides the sufficient condition for a localextremum.The analysis so far gives that a local maximum of ψ ( p ) ei-ther has: (i) ¯ ψ ( p , p ′ ) = 0 in any direction of perturbation,or (ii) ¯ ψ ( p , p ′ ) < for all feasible perturbations. It is notpossible to find a local maximum such that the linear varia-tion vanishes in some directions of perturbation and negativein others. Notice that case (i) only provides a necessary con-dition for a local extremum and we need to further determineif it is a maximum or a minimum.To proceed, we continue to analyze the local extrema with ¯ ψ ( p , p ′ ) = 0 from Eq. (21) which, for x ∗ i > , can be rewrit-ten as n X j =1 p j k j = 2 n X j =1 p j k j x ∗ j . (22)The left-hand side equals z whereas the right side equals z Θ ∗ , implying Θ ∗ = 1 / . Since, at equilibrium, we have x ∗ i = λk i Θ ∗ λk i Θ ∗ = λk i λk i , (23)from the definition of Θ ∗ , we obtain the following relation fora local extremum: z Θ ∗ = n X i =1 p i k i λk i λk i = z . (24)Together with the probability and the average degree con-straints, a local extremum with ¯ ψ ( p , p ′ ) = 0 can be foundin the set P o where any p ∈ P o satisfies n X i =1 p i λk i λk i = z , n X i =1 p i k i = z, n X i =1 p i = 1 , p i ∈ [0 , ∀ i ∈ { , · · · , n } . (25)The spreading prevalence for p ∈ P o can be directly obtainedfrom the definition of P o , without solving the HMF equations.In particular, subtracting the second equation in Eq. (25) bythe first equation on both sides, we have n X i =1 p i k i λk i = 2 λ n X i =1 p i x i = 2 λ ψ ( p ) = z , (26)which implies ψ ( p ) = λz/ for p ∈ P o . That is, for any p ∈ P o , the resulting spreading prevalence is the same.For p ∈ P o , conversely we have ¯ ψ ( p , p ′ ) = 0 for all feasi-ble directions. To see this, consider the definition of Θ ∗ , z Θ ∗ = n X i =1 p i k i λk i Θ ∗ λk i Θ ∗ . (27)If the right-hand side is viewed as a function of Θ ∗ , then itincreases with Θ ∗ . For Θ ∗ = 0 , the right-hand side equalszero and for Θ ∗ → ∞ it converges to z . Consequently, fora fixed p , there is a unique Θ ∗ such that the right-hand sideequals z/ . Since p ∈ P o , from the first equation in Eq. (25),we have Θ ∗ = 1 / and then Eq. (21) holds. Similarly, for p / ∈ P o , we have Θ ∗ = 1 / . The conclusion is that for p ∈ P , ¯ ψ ( p , p ′ ) = 0 holds for all feasible directions if andonly if p ∈ P o . C. Necessary condition for the global optimal solution
Suppose P o is nonempty, the question is as follows: Arethe degree distributions local maxima or a global maximum?As the set P o is defined through simple linear equations, wecan prove that any p ∈ P o is indeed a global maximum via al-gebraic manipulations. Concretely, in the following, we provethat if p / ∈ P o , then ψ ( p ) < λz/ . When P o is empty, weshow that the support of the optimal degree distribution hasno more than two distinct nodal degrees.For any p / ∈ P o , this is trivially true if Θ ∗ = 0 and weassume Θ ∗ > . Suppose there exists p / ∈ P o but ψ ( p ) ≥ λz/ , then from the definition of ψ ( p ) , we have λ Θ ∗ ψ ( p ) = n X i =1 p i k i λk i Θ ∗ ≥ z ∗ . (28)Subtracting P ni =1 p i k i = z from the inequality on both sides,we have n X i =1 p i k i λk i Θ ∗ λk i Θ ∗ = z Θ ∗ ≤ z − z ∗ . (29) The inequality implies (2Θ ∗ − ≤ . An equality holdsonly when Θ ∗ = 1 / , but this contradicts with Θ ∗ = 1 / for p / ∈ P o from the discussions below Eq. (27).The analysis so far reveals that, when P o is nonempty, any p is a global maximum if and only if it belongs to P o . It re-mains to address the following issues. (i) For which values of λ is the set P o nonempty? (ii) If P o is empty, how do we findthe local maxima with ¯ ψ ( p , p ′ ) < for all feasible perturba-tions. We will solve (ii) partly for the rest of this section, andprovide full answers to (i) and (ii) in the next section.Suppose P o is empty. Consider any p ∈ P and definethe support of p as supp( p ) = { k i : p i > } . Suppose supp( p ) has more than two distinct nodal degrees, we canpick any m ≥ points { k i , k i , · · · k i m } ⊂ supp( p ) fromthe support of p and consider a perturbation ¯ p whose elementsare nonzero only at these points. For any ¯ p which is nonzeroonly on the support of p , we can always choose α sufficientlysmall such that p + α ¯ p ∈ [0 1] n , p − α ¯ p ∈ [0 1] n . (30)The perturbations α ¯ p and − α ¯ p are thus both feasible for suf-ficiently small α . As P o is empty, there always exists ¯ p suchthat Z ( m ) ¯ p = 0 . From Eq. (17), we have ¯ ψ ( p , p + α ¯ p ) = − ¯ ψ ( p , p − α ¯ p ) . (31)This indicates that if P o is empty, then any p whose supporthas more than two distinct degrees cannot be a local maximumand the optimal p o must concentrate on no more than twodistinct nodal degrees. IV. FINDING THE OPTIMAL DEGREE DISTRIBUTIONS
The results in Sec. III indicate that, to find the optimal dis-tributions, it is only necessary to determine whether set P o isnonempty. If it is empty, the task is to search through all de-gree distributions whose support consists of one or two nodaldegrees. In fact, in the latter case, the HMF equation can besolved analytically to yield the optimal degree distributions. A. Conditions for P o to be nonempty As P o is a closed convex set, by the Krein-Milman the-orem, it is the convex hull of all its extremum points (i.e., p ∈ P o that does not lie in the open line segment joiningany two other points in P o ). To check if P o is nonempty isequivalent to examining if all its extremum points exist. Inthe following, we first show that the support of the extremumpoints of P o has no more than three distinct nodal degrees.In this case, the value of p is uniquely determined by choiceof the support. As a result, we can solve p in terms of thesupport and λ explicitly. With a fixed chosen support and the λ value so determined, the corresponding p is physical for p ∈ [0 , n . By checking all the points that are supported onno more than three degrees, we can derive the condition for λ under which P o is nonempty.Suppose there exists p ∈ P o whose support has more thanthree degrees. Pick any m ≥ points { k i , k i , · · · k i m } ⊂ supp( p ) and consider a perturbation ¯ p whose elements arenonzero only on these points. Define Y ( m ) = · · · k i k i · · · k i m λk i λk i λk i λk i · · · λk im λk im . (32)A feasible direction of perturbation ¯ p , which keeps p ± α ¯ p staying inside P o for sufficiently small values of α , mustsatisfy the condition Y ( m ) ¯ p = 0 . The nullity of Y ( m ) is nullity( Y ( m ) ) = m − . Thus, for m > , the space offeasible perturbations is nonempty. Moreover, we can alwayschoose α > and α > such that the support of p + α ¯ p and p − α ¯ p has m − distinct nodal degrees. In this way, p lies on the open line segment that joins p + α ¯ p and p − α ¯ p .This means that, if the support of p ∈ P o has more than threedistinct nodal degrees, it will not be an extremum point of P o .To determine if P o is nonempty, it thus suffices to checkif there exists p ∈ P o whose support has no more than thanthree distinct nodal degrees. Consider any k i > k i > k i ,the values of p i , p i and p i are uniquely determined byEq. (25), which are p i = − ( k i λ + 2) g ( k i , k i )8 λ ( k i − k i )( k i − k i ) ,p i = + ( k i λ + 2) g ( k i , k i )8 λ ( k i − k i )( k i − k i ) ,p i = − ( k i λ + 2) g ( k i , k i )8 λ ( k i − k i )( k i − k i ) , (33)where g ( k a , k b ) = (cid:0) λ z − λ (cid:1) k a k b + 2 λz ( k a + k b ) − z. (34)The degree distribution is physically meaningful insofar as p i , p i , p i ∈ [0 1] . Since p i + p i + p i = 1 , it is suffi-cient to guarantee p i , p i and p i to be nonnegative, i.e., toguarantee g ( k i , k i ) ≤ , g ( k i , k i ) ≤ , g ( k i , k i ) ≥ . (35)In Appendix D, we analyze the three inequalities in detail.Here we summarize the procedure and results. We study un-der what conditions the three inequalities in (35) hold con-secutively. Particularly, we first derive the condition for theexistence of ( k i , k i ) such that g ( k i , k i ) ≥ holds. Then,under this condition, we check if there exists k i such that theother two inequalities in (35) hold. Consider the inequality g ( k i , k i ) ≥ . The possible values of the two nodal degreesare k i ∈ { k , k , · · · , z + } and k i ∈ { z − , · · · , k n − , k n } ,where z + = min i { k i ≥ z } and z − = max i { k i ≤ z } . As g ( k a , k b ) is quadratic in λ , we can show that g ( k a , k b ) ≥ if λ ≥ λ ( k a ,k b ) but g ( k a , k b ) < otherwise, where λ ( k a ,k b ) = 2 z − k a − k b + s(cid:18) k a + 1 k b − z (cid:19) + 4 k a k b . (36) As λ ( k a ,k b ) is a decreasing function of k a for k a ≥ z + and anincreasing function of k b for k b ≤ z − , we can show that thereexists ( k i , k i ) such that g ( k i , k i ) ≥ holds insofar as λ ≥ λ , where λ = λ ( k ,k n ) . Furthermore, when this conditionholds, we can show that there exists k i such that the other twoinequalities in (35) hold if and only if λ ≤ λ = λ ( z + ,z − ) .Overall, the values of λ are divided by λ , λ , and λ intofour regions, where λ is defined in Eq. (5). The four regionsare described as follows.(i) For λ ≤ λ , the optimization problem is trivial, i.e., nodegree distribution can trigger an epidemic outbreak.(ii) For λ < λ < λ , set P o is empty, thus the global max-imum can only be found among all p supported on one or twonodal degrees.(iii) For λ ≤ λ ≤ λ , set P o is nonempty and any p ∈ P o will lead to equal spreading prevalence λz/ . In Appendix D,we show that for λ = λ , set P o consists of a unique degreedistribution supported on { k , k n } , whereas for λ = λ , set P o has a unique degree distribution supported on { z + , z − } .For λ < λ < λ , there are infinitely many global maximathat constitute a plateau of equal spreading prevalence.(iv) For λ > λ , set P o again becomes empty, and theglobal maxima can only be supported on one or two nodaldegrees. B. Analytic solutions of HMF equations on one or two degrees
Having determined the conditions under which P o isnonempty, we are now in a position to find the optimal de-gree distributions that are supported on one or two degrees.In this case, the HMF equations consist of only one or twodifferent equations so the equilibrium solution can be solvedexplicitly. We can then directly optimize the solution to obtainthe optimal degree distribution on one or two nodal degrees.Consider the situation where p is supported on one or twodifferent nodal degrees. Let k ≥ k i ≥ z + and z − ≥ k i ≥ k n be any two nodal degrees from k so that p i and p i areuniquely determined by p i + p i = 1 , p i k i + p i k i = z, (37)which leads to the solutions of p i and p i in terms of k i , k i ,and z as p i = z − k i k i − k i , p i = k i − zk i − k i . (38)When z is an integer and either k i or k i equals z , it reducesto the case where p is supported on one nodal degree. Withthe values of p i and p i , the HFM equation can be solved an-alytically (Appendix E). After some algebraic manipulations,we obtain the spreading prevalence as ψ ( p ) =1 − u − λz (cid:0) u + v (cid:1) + uλz p λz ( λz − u ) + 4 v . (39)where u = 12 (cid:18) zk i + zk i (cid:19) , v = 12 (cid:18) zk i − zk i (cid:19) . (40)The degrees k i and k i are then uniquely determined by thevalues of u and v .We can now carry out optimization among all degree distri-butions that are supported on one or two nodal degrees. Thegoal is to find the optimal degree values k i and k i such that ψ ( p ) given by Eq. (39) is maximized. Our approach is totreat k i and k i as continuous variables to obtain the maximaof ψ ( p ) , which can finally be used to find the actual optimalvalues of k i and k i as integers.From Eq. (38), we see that p i and p i are uniquely deter-mined by the choice of k i and k i which, in turn, are uniquelydetermined by the values of u and v defined in Eq. (40). Theequivalent problem is to optimize ψ ( p ) by u and v . It is con-venient to rewrite ψ ( p ) as ψ ( u, v ) . Taking the partial deriva-tives of ψ ( u, v ) , we obtain ∂ψ ( u, v ) ∂u = (cid:18) λz p λz ( λz − u ) + 4 v − (cid:19) × − u p λz ( λz − u ) + 4 v ! , (41) ∂ψ ( u, v ) ∂v = 2 vλz u p λz ( λz − u ) + 4 v − ! . (42)The two partial derivatives vanish simultaneously only for u = p λz ( λz − u ) + 4 v , (43)which defines a curve on the u - v plane where every pointon it is a critical point of ψ ( u, v ) . Substituting Eq. (43) intoEq. (39), we obtain the spreading prevalence along the curveas ψ ( p ) = λz , (44)which is exactly the spreading prevalence for those p ∈ P o ,given that P o is nonempty.Substituting the definition of u and v in Eq. (40) intoEq. (43), we can express the curve in terms of k i and k i as g ( k i , k i ) = 0 , where g ( k a , k b ) = ( λ z − λ ) k a k b + 2 λz ( k a + k b ) − z. (45)This function is also exactly the same as Eq. (34), the onethat emerges when we analyze the extremum points of P o .Not all points ( k a , k b ) along the optimal curve in Eq. (43) arephysically meaningful. Especially, for a point on the k a - k b plane to be meaningful, it must be an integer point that lies inthe region, R = { ( k a , k b ) : k ≥ k a ≥ z + , z − ≥ k b ≥ k n } . (46)From the discussions below Eq. (34), the curve g ( k a , k b ) = 0 passes an integer point ( k a , k b ) when λ = λ ( k a ,k b ) , where λ ( k a ,k b ) is defined in Eq. (36). When this happens, the degreedistribution supported on { k a , k b } belongs to set P o . In fact,if we let ( k i , k i ) = ( k a , k b ) and substitute g ( k i , k i ) = 0 into Eq. (33), we then have p i = 0 and p i = z − k i k i − k i , p i = k i − zk i − k i . (47)This recovers exactly the same degree distribution defined inEq. (38). For λ < λ or λ > λ , set P o is empty, and nointeger point in region R can lie on the curve g ( k a , k b ) = 0 . Inthis case, it is necessary to further analyze the optimal degreedistribution.For convenience, we write ψ ( p ) as ψ ( k i , k i ) and have ∂ψ ( k i , k i ) ∂k i = − z k i (cid:18) ∂ψ ( u, v ) ∂u + ∂ψ ( u, v ) ∂v (cid:19) . (48)Substituting Eqs. (41) and (42) into Eq. (48), we get ∂ψ ( k i , k i ) ∂k i = − u p λz ( λz − u ) + 4 v ! × z k i (cid:18) vλz + 1 − λz p λz ( λz − u ) + 4 v (cid:19) . (49)Since u − v = z/k i > , we have p λz ( λz − u ) + 4 v > p λz ( λz + 4 v ) + 4 v > v + λz. (50)The last line in Eq. (49) is, thus, negative. For p λz ( λz − u ) + 4 v − u > , (51) ψ ( k i , k i ) is a decreasing function of k i ; otherwise it is anincreasing function of k i . Similarly, the partial derivative of ψ ( k i , k i ) with respect to k i is ∂ψ ( k i , k i ) ∂k i = − u p λz ( λz − u ) + 4 v ! × z k i (cid:18) − vλz + 1 − λz p λz ( λz − u ) + 4 v (cid:19) , (52)where the term in the last line is positive. Thus, if Eq. (51)holds, ψ ( k i , k i ) is an increasing function of k i , otherwiseit is a decreasing function of k i .Recall that g ( k a , k b ) in Eq. (34) is equivalent to the relationin Eq. (43). The inequality in Eq. (51) can then be written interms of k i and k i as g ( k i , k i ) > . (53)From the discussions in Appendix D, for any ( k a , k b ) ∈ R ,we have g ( k a , k b ) < if λ < λ and g ( k a , k b ) > if λ > λ .These results lead to the optimal degree distributions in eachof the parameter regions of λ .For λ < λ < λ , we have g ( k a , k b ) < for any ( k a , k b ) ∈ R . Consequently, ψ ( k i , k i ) is an increasing func-tion of k i and a decreasing function of k i . In this case,the optimal degree distribution is supported on k a = k and k b = k n . Moreover, the spreading prevalence of the optimaldegree distribution is strictly less than λz/ .For λ ≤ λ ≤ λ , the degree distribution p is a globalmaximum if and only if p ∈ P o . Since P o is a connected set,all the global maxima constitute a plateau of degree distribu-tions with equal spreading prevalence. For λ = λ , the set P o consists of a unique degree distribution, which is exactlythe optimal one for λ < λ . For λ = λ , the set P o also hasone unique degree distribution, and we will see that it is theoptimal one for λ > λ .For λ > λ , we have g ( k a , k b ) > for any ( k a , k b ) ∈ R .As a result, ψ ( k i , k i ) is a decreasing function of k i and anincreasing function of k i . In this case, the optimal degreedistribution is supported on k i = z + and k i = z − . V. CHARACTERISTICS OF OPTIMAL DEGREEDISTRIBUTIONS
For relatively low infection rates ( λ < λ ≤ λ ), the op-timal degree distribution is supported on the maximal andminima possible degrees { k , k n } . For high infection rates( λ ≥ λ ), the optimal degree distribution is supported on thetwo nodal degrees { z + , z − } that are nearest to the averagedegree z . Therefore, we need to study how the support of theoptimal degree distributions behaves for intermediate infec-tion rates in the range [ λ , λ ] . Let P e ⊂ P o be the set of allextremum points of P o , where P e is a finite set. As P o is theconvex hull of all its extremum points, for any p ∈ P o , it is aconvex combination of the extremum points, p = X p e ∈P e c ( p e ) p e , (54)where c ( p e ) ≥ and P p e ∈P e c ( p e ) = 1 . The broadestsupport (i.e., the support with the largest number of distinctnodal degrees) of p ∈ P o thus is [ p e ∈P e supp ( p e ) . (55)Any degree distribution with c ( p e ) > for all p e ∈ P e willhave the broadest possible support.Consider the case where λ is slightly above λ and ( k , k n ) is a unique point such that g ( k , k n ) > . From the discus-sions at the end of Appendix D, we have that, by choosing k i = k , k i = k n , and k i to be any allowed degree with k > k i > k , the triple ( k i , k i , k i ) will define a physicaldegree distribution from Eq. (33). As the middle degree k i isarbitrary, the broadest support in this case consists of all theallowed degrees in k , i.e., the cardinality of the broadest sup-port increases abruptly from 2 to n at λ = λ . Similarly, itcan be seen that, when λ is slightly below λ and ( z + , z − ) isthe unique point such that g ( z + , g − ) < , the broadest sup-port also consists of all the possible nodal degrees. Figure 1 FIG. 1. Normalized cardinality of the broadest support for p ∈ P o versus λ . The vertical gray dashed lines mark the locations of λ , λ ,and λ that divide the values of λ into different regions. For λ ≤ λ or λ ≥ λ , the normalized cardinality is /n , whereas it is one for λ < λ < λ . For λ < λ < λ , the cardinality of the broad-est possible support is obtained by testing all the extremum pointsnumerically. The values of other parameters are k = 30 , k n = 1 ,and z = 15 . . The values of λ i for i ∈ { – } are λ ≈ . , λ ≈ . and λ ≈ . . The support of the degree dis-tribution can take on any integer value between k and k n , i.e., k = [30 , , · · · , T and n = 30 . shows the normalized cardinality of the broadest possible sup-port versus λ . We see that, for λ < λ < λ , the broadestsupport indeed consists of all the distinct degrees allowed in k , indicating that, except for relatively low or high values of λ , the support of the optimal degree distribution can be quitebroad.In general, the degree heterogeneity of a network, definedas H = h k i / h k i , can have significant impacts on thespreading dynamics. A natural question is, what is the de-gree heterogeneity of the optimal degree distribution? Sincethe average degree is fixed ( h k i = z ), the degree heterogene-ity determines the outbreak threshold. For sufficiently smallvalues of λ where there is a unique network that can trigger anepidemic outbreak, the optimal network structure is one withthe largest degree heterogeneity.Consider the general problem of finding maxima and min-ima of H among all degree distributions. The extrema of H can be found by maximizing or minimizing the secondmoment h k i of the degree distribution. The Bhatia-Davisinequality stipulates that the second moment of p is max-imized when it is concentrated at the endpoints k and k n .To minimize the second moment, we note that the definition h k i = P ni =1 p i k i has a similar form to Eq. (17) with χ i re-placed by k i . Following the reasoning in Sec. III B, we seethat the minimum of H is supported on two nodal degrees.Through a direct comparison of all distributions supported ontwo degrees, we find that H is minimized when p concentrateson { z + , z − } . We see that the optimal degree distributions for FIG. 2. Bounds of degree heterogeneity of the optimal degree distri-butions versus the infection rate. The vertical gray dashed lines markthe locations of λ , λ , and λ that divide λ into different regions.The blue solid trace represents the bounds of the degree heterogene-ity H . For λ ≤ λ or λ ≥ λ , the lower and upper bounds coincide.For λ < λ < λ , the degree heterogeneity can take on any valuein the shaded region. Other parameter values are k = 30 , k n = 1 and z = 15 . . The values of λ i for i ∈ { – } and n are the same asthose in Fig. 1. λ ≤ λ and λ ≥ λ are exactly the ones that maximize andminimize the degree heterogeneity, respectively.For a fixed λ value in the intermediate region ( λ < λ <λ ), the values of H for different degree distributions in P o are not necessarily identical. From Eq. (54), we see that, if p is a convex combination of the extremum points, its secondmoment can be obtained by the same convex combination ofthe second moment of the extremum points. Consequently,the degree heterogeneity of p ∈ P o is bounded by that of theextremum points. Figure 2 shows the bounds of the degreeheterogeneity H of the optimal degree distributions versus λ .The general phenomenon is that the optimal network is moreheterogeneous for small infection rates but less so for largerates, as the upper and lower bound of H decreases with λ .However, the degree heterogeneity does not decrease with λ in a strict sense but only trendwise. In fact, if we draw a linesegment joining the two degree distributions that reach the up-per and lower bounds, then H varies continuously on this linesegment, i.e., the degree heterogeneity can take on any valuebetween the lower and upper bounds.Our analysis of the characteristics of the optimal degree dis-tributions reveals a phenomenon: The existence of a particularvalue of the infection rate for which every degree distributionis optimal. From the definition of P o , any p ∈ P o must sat-isfy the first equation in Eq. (25), whose left-hand side is an in-creasing function of λ that converges to zero or z for λ → or λ → ∞ , respectively. As a result, for any p ∈ P , there alwaysexists a unique λ value such that p ∈ P o . Only two degreedistributions are optimal under multiple values of λ , which arethe two supported on either { k , k n } or { z + , z − } , as they are ( + )/2 FIG. 3. Spreading prevalence divided by λz versus λ for three dif-ferent degree distributions. The values of the spreading prevalenceare obtained by solving the HMF equations numerically. The ver-tical gray dashed lines mark the locations of λ , ( λ + λ ) / , and λ where the three degree distributions are optimal. The horizontalblack dashed line correspond to ψ ( p ) /λz = 1 / . Other parametervalues are k = 30 , k n = 1 , and z = 15 . . The values of λ i for i ∈ { – } and n are the same as those in Fig. 1. optimal when P o is empty. In Sec. III B, we have shown thatfor any p / ∈ P o , its spreading prevalence is strictly less than λz/ . This suggests the following phenomenon: For any de-gree distributions, its spreading prevalence as a function of λ will touch the line ψ ( p ) = λz/ only at one value of λ andunder this value of λ the degree distribution is among the opti-mal degree distributions. For all other values of λ , its spread-ing prevalence will strictly be below the line ψ ( p ) = λz/ .To illustrate the phenomenon, we consider three degree dis-tributions that are optimal at λ = λ , λ = ( λ + λ ) / , and λ = λ , respectively. For λ = λ or λ = λ , the optimaldegree distribution is unique. For λ = ( λ + λ ) / , we ran-domly pick a degree distribution from P o by a uniformly ran-dom convex combination of the extremum points. We plot ψ ( p ) /λz versus λ for three degree distributions as shown inFig. 3. It can be seen that the value of ψ ( p ) /λz reaches 1/4 atthe predicted value of λ and is below 1/4 for any other valuesof λ . VI. DISCUSSION
Given a dynamical process of interest, identifying the ex-tremum network provides deeper insights into the interplaybetween network structure and dynamics. From the perspec-tive of applications, searching for a global dynamics-specificoptimal network can be valuable in areas such as informationdiffusion, transportation, and behavior promotion. The issue,however, belongs to the category of dynamics-based inverseproblems that are generally challenging and extremely diffi-cult to solve. We have taken an initial step in this direction.0Specifically, by limiting the study to SIS type of spreadingdynamics and imposing the annealed approximation, we haveobtained analytic solutions to the inverse problem. Our solu-tions unveil a phenomenon with implications: A fundamentalcharacteristic of the optimal network, its degree heterogene-ity, depends on the infection rate. In particular, strong degreeheterogeneity facilitates the spreading but only for small in-fection rates. For relatively large infection rates, the optimalstructure tends to choose the networks that are less heteroge-neous. This means that, when designing an optimal network,e.g., for information spreading, the ease with which informa-tion can diffuse among the nodes must be taken into account.Our analysis has also revealed the existence of a particularvalue of the infection rate for which every degree distributionis globally optimal.The annealed approximation that serves the base of ouranalysis is applicable to networks that are describable by theuncorrelated configuration model. It remains to be an openproblem to find the optimal quenched networks for SIS dy-namics. In Ref. [16], the authors introduced a technique tobridge the annealed and quenched limit of the SIS model. Thetechnique can provide a starting point to extend our analyticapproach to quenched networks. The variational analysis inthe current paper can be extended to SIS type dynamics onquenched weighted networks to derive a necessary conditionfor local optimum. In the variational calculus, we have to per-form a network structural perturbation to the mean-field equa-tion; therefore, we emphasize a necessary element that makesthe variational calculations viable: The spreading prevalenceis a continuous function of the perturbations, at least locallyaround the network being perturbed. The variational analysiswill result in a necessary condition for local optimum. How-ever, it is not clear yet what we can derive from the neces-sary condition without annealed approximations. To general-ize the theory to settings under less stringent simplificationsis at present an open topic worth investigating. Another as-sumption in the present paper is that only the number of edgesis held fixed, and it is useful to study the optimal networksunder more realistic restrictions. Moreover, it is of general in-terest to seek optimal solutions of network structures for dif-ferent types of dynamical processes. Our paper represents astep forward in this direction.
ACKNOWLEDGMENTS
L.P. would like to acknowledge support from the Na-tional Natural Science Foundation of China under GrantNo. 62006122. W.W. would like to acknowledge supportfrom the National Natural Science Foundation of China underGrant No. 61903266, Sichuan Science and Technology Pro-gram Grant No. 20YYJC4001, China Postdoctoral ScienceSpecial Foundation Grant No. 2019T120829, and the Fun-damental Research Funds for the central Universities GrantNo. YJ201830. L.T. would like to acknowledge support fromThe Major Program of the National Natural Science Founda-tion of China under Grant No. 71690242 and the National KeyResearch and Development Program of China under Grant No. 2020YFA0608601. Y.-C.L. would like to acknowledgesupport from the Vannevar Bush Faculty Fellowship programsponsored by the Basic Research Office of the Assistant Sec-retary of Defense for Research and Engineering and fundedby the Office of Naval Research through Grant No. N00014-16-1-2828.
Appendix A: Proof that x ∗ is a continuously differentiablefunction of p above the outbreak threshold Setting the right-hand side of Eq. (1) to zero for an equilib-rium, we get x ∗ i = λk i Θ ∗ λk i Θ ∗ , (A1)where Θ ∗ is obtained by substituting x ( t ) = x ∗ into Eq. (2).Define a function f ( p , x ) : ( p , x ) → R n as f i ( p , x ) = x i − λk i P nl =1 p l k l x l z + λk i P nl =1 p l k l x l , (A2)we have f ( p , x ∗ ) = 0 . Note that f ( p , x ) is a continuouslydifferentiable function of p and x . We now show that, fromthe relation f ( p , x ∗ ) = 0 , the stable equilibrium point x ∗ canbe written as a continuously differentiable function of p for λ > z/ h k i by applying the implicit function theorem.The derivative of f i ( p , x ) with respect to x j is ∂f i ( p , x ) ∂x j = δ i,j − λzk i k j p j ( z + λk i P nl =1 p l k l x l ) , (A3)where δ i,j is the Kronecker δ . The Jacobian matrix of f ( p , x ) to x can be written as I − rs T , where I is the n × n identitymatrix and r and s are n × vectors with elements, r i = λzk i ( z + λk i P nl =1 p l k l x l ) , s i = k i p i . (A4)Let b be an eigenvector of the matrix rs T with eigenvalue ω .From the eigenvalue equation, we have n X j =1 r i s j b j = ωb i . (A5)Multiplying both sides by s i and summing over i , we have n X i =1 r i s i n X j =1 s j b j = ω n X i =1 s i b i . (A6)As a result, the only possible eigenvalues of matrix rs T are ω = 0 or ω = P ni =1 r i s i .For λ > z/ h k i , all elements of x ∗ are positive. At x = x ∗ ,1we have n X i =1 r i s i = n X i =1 λzp i k i ( z + λk i P nl =1 p l k l x ∗ l ) = 1 z n X i =1 λp i k i (1 − x ∗ i ) = 1 z Θ ∗ n X i =1 p i k i x ∗ i (1 − x ∗ i )=1 − z Θ ∗ n X i =1 p i k i ( x ∗ i ) < , (A7)where the second and third equalities can be verified by sub-stituting them into Eq. (A1) and x ∗ i = λk i (1 − x ∗ i )Θ ∗ , respec-tively.Taken together, the eigenvalues of the matrix rs T are lessthan one for λ > z/ h k i , so the eigenvalues of the Jacobianmatrix I − rs T are less than zero, which further implies thatthe Jacobian matrix is invertible. By the implicit function the-orem, x ∗ is a continuously differentiable function of p . Appendix B: Eigenvalues of the Jacobian matrix
Denote the right side of Eq. (1) by h i = − x i ( t ) + λk i [1 − x i ( t )] Θ . (B1)The Jacobian matrix for h = ( h , · · · , h n ) T at x = x ∗ isexactly J : ∇ h = J . As x ∗ is the unique global stable equi-librium point [15], the eigenvalues of J must have negativereal parts. Appendix C: Detailed derivation of χ Define two vectors µ and ν of length n whose elements are µ i = λz k i (1 − x ∗ i ) , ν i = k i p i . (C1)Further, define a n × n diagonal matrix D with the elements D ii = − − λk i Θ ∗ . (C2)By the Sherman-Morrison formula, we have J − = (cid:0) D + µ · ν T (cid:1) − = D − − D − · µ · ν T · D − ν T · D − µ . (C3)Substituting Eq. (14) into Eq. (16), we get Eq. (17) with χ i given by χ i = x ∗ i − k i x ∗ i p T J − µ. (C4)Inserting Eq. (C3) into Eq. (C4) leads to χ i = x ∗ i + λk i x i P nj =1 p j k j (1 − x ∗ j ) (1 + λk j Θ ∗ ) − z − λ P nj =1 p j k j (1 − x ∗ j ) (1 + λk j Θ ∗ ) − . (C5) At the equilibrium point, we have − x ∗ i + λk i (1 − x ∗ i ) Θ ∗ = 0 , (C6)which leads to (1 + λk j Θ ∗ ) − = (1 − x ∗ i ) . (C7)Substituting the above two equations into Eq. (C5), we obtain χ i = x ∗ i λk i Θ ∗ P nj =1 p j k j (1 − x ∗ j ) P nj =1 p j k j ( x ∗ j ) ! , (C8)which is Eq. (18). Appendix D: Conditions for P o to be nonempty We test the validity of the three inequalities in (35) in asequential manner: First we study the condition for λ whenthere exist k i and k i such that g ( k i , k i ) ≥ holds, wethen test under the derived condition if there exists k i suchthat the other two inequalities hold.As a preparation, we prove a result that will be used re-peatedly in the rest of this appendix. In particular, we showthat for P o to be nonempty, it is necessary to have λz ≤ from Eq. (26). Note that Eq. (26) is the average of thefunction f ( k a ) = 2 k a / (2 + λk a ) under the degree distribu-tion p . This function has a negative second order derivative f ′′ ( k a ) = − λ/ (2 + λk a ) , so f ( k a ) is concave. By Jensen’sinequality, we have n X i =1 p i k i λk i ≤ z λz . (D1)Since the left side equals z/ from Eq. (26), it is necessary tohave z/ (2 + λz ) ≥ z/ , which implies λz ≤ . The equalityin Eq. (D1) holds only when z is an integer and is one of theallowed degrees in k and, in addition, p concentrates on z . Inthis case we have λz = 2 , so P o has a unique element p thatconcentrates on z .We consider the case of λz < . The analysis be-gins with the setting of the existence of ( k i , k i ) such that g ( k i , k i ) ≥ holds. Defining z + = min i { k i ≥ z } and z − = max i { k i ≤ z } , we have k i ∈ { k , k , · · · , z + } and k i ∈ { z − , · · · , k n − , k n } . The function g ( k a , k b ) isquadratic in λ , and the equation g ( k a , k b ) = 0 has two roots:one positive and one negative. The positive one is λ ( k a ,k b ) = 2 z − k a − k b + s(cid:18) k a + 1 k b − z (cid:19) + 4 k a k b . (D2)As a result, for < λ < λ ( k a ,k b ) , we have g ( k a , k b ) < ,whereas g ( k a , k b ) ≥ for λ ≥ λ ( k a ,k b ) . If λ ( k a ,k b ) is re-garded as a function of k a and k b , through the derivatives, wehave that λ ( k a ,k b ) is a decreasing function of k a for k a > z and an increasing function of k b for k b < z . Consequently,the value of λ ( k a ,k b ) reaches its minimum at ( k , k n ) . There2exists at least one ( k i , k i ) such that g ( k i , k i ) ≥ insofaras λ ≥ λ ( k ,k n ) .Having determined the condition under which there exists ( k i , k i ) such that g ( k i , k i ) ≥ holds, we can obtain theconditions under which there exists k i such that g ( k i , k i ) ≤ and g ( k i , k i ) ≤ . When the curve g ( k a , k b ) = 0 passesan integer point that can be chosen as ( k i , k i ) , we have g ( k i , k i ) = 0 . From Eq. (33), we have p i = 0 and p i = z − k i k i − k i , p i = k i − zk i − k i . (D3)We see that p i and p i are independent of the choice of k i and p is supported on one or two nodal degrees.Now consider the case of g ( k i , k i ) > . For fixed ( k i , k i ) , the inequalities g ( k i , k i ) ≤ and g ( k i , k i ) ≤ can be rearranged as (cid:2)(cid:0) λ − λ z (cid:1) k i − λz (cid:3) k i ≥ λzk i − z, (cid:2)(cid:0) λ − λ z (cid:1) k i − λz (cid:3) k i ≥ λzk i − z. (D4)As k i ≥ z and λz < , we have (cid:0) λ − λ z (cid:1) k i − λz > . (D5)From g ( k i , k i ) ≥ we have (cid:0)(cid:0) λ − λ z (cid:1) k i − λz (cid:1) k i ≤ λzk i − z. (D6)Because λz < and k i ≤ z , the right-hand side of the aboveinequality is negative. We, thus, have (cid:0) λ − λ z (cid:1) k i − λz < . (D7)With the above results, Eq. (D4) implies that there exist feasi-ble values of k i insofar as λzk i − z (4 λ − λ z ) k i − λz ≤ λzk i − z (4 λ − λ z ) k i − λz (D8)and there is at least one integer between the two sides of theinequality.Defining a function of λ and k a as f ( λ, k a ) = 2 λzk a − z (4 λ − λ z ) k a − λz , (D9)we have that the left and right sides of Eq. (D8) are equal to f ( λ, k i ) and f ( λ, k i ) , respectively. The derivative of f ( k a ) with respect to k a is ∂f ( λ, k a ) ∂k a = 8 λz (2 − λz )((4 λ − λ z ) k a − λz ) . (D10)Consequently, f ( λ, k a ) is an increasing function of k a for λz < and the function is non-singular, so f ( λ, k i ) isbounded from above as f ( λ, k i ) < lim k a →∞ f ( k a ) = 2 z − λz < z, (D11) whereas f ( λ, k i ) is bounded from below as f ( λ, k i ) > lim k a → f ( k a ) = 2 λ > z. (D12)We thus have that the inequality f ( λ, k i ) < f ( λ, k i ) holdsfor λz < . It remains to determine if there is an integer be-tween f ( λ, k i ) and f ( λ, k i ) . In this regard, if z is an integerand is one of the degrees allowed, the situation is relativelysimple, and we pick k i = z .We analyze the case where z is not an integer. Note thatthe left-hand side of Eq. (D8) is strictly less than the right-hand side and f ( λ, k a ) is an increasing function of k a . Thegap between the two sides of Eq. (D8) is then maximized for ( k i , k i ) = ( z + , z − ) . Suppose there are no integer pointsbetween f ( λ, z + ) and f ( λ, z − ) . It implies that there are nointeger points between f ( λ, k i ) and f ( λ, k i ) for any otherchoice of ( k i , k i ) . For λ = λ ( z + ,z − ) , the curve g ( k a , k b ) =0 passes the point ( k a , k b ) = ( z + , z − ) and set P o is nonemptybased on Eq. (D3), as we can take ( k i , k i ) = ( z + , z − ) . Inthe next, we show that if λ > λ ( z + ,z − ) , set P o will be emptyas there are no integer points between f ( λ, z + ) and f ( λ, z − ) .However, for λ < λ ( z + ,z − ) , P o is guaranteed to be nonempty.To show that P o is empty for λ > λ ( z + ,z − ) , we note thatthe derivative of f ( λ, k a ) with respect to λ is ∂f ( λ, k a ) ∂λ = 2 z ( λk a − + 16 z ( k a − z )[(4 λ − λ z ) k a − λz ] . (D13)For k a ≥ z and λz < , the derivative is positive, so f ( λ, k a ) is an increasing function of λ . Now we show that if (cid:0) λ − λ z (cid:1) k a − λz < , (D14)then f ( λ, k a ) is a decreasing function of λ . Note that k i sat-isfies the above inequality for g ( k i , k i ) > [c.f., the discus-sions above Eq. (D8)]. Taking the derivative with respect to k a for the numerator of the right-hand side of Eq. (D13), weget λz ( λk a −
2) + 16 z > z − λz > , (D15)where the second inequality is the result of applying λz < . The numerator on the right side of Eq. (D13) itself is anincreasing function of k a . In addition, Eq. (D14) implies k a < z − λz . (D16)When k a equals the right side of this inequality, the numeratorof the right side of Eq. (D13) becomes, λz ( λz − − λz ) < , (D17)so f ( λ, k a ) is a decreasing function of λ when Eq. (D14)holds. For λ = λ ( z + ,z − ) , we have f ( λ, z + ) = z − and f ( λ, z − ) = z + . For λ > λ ( z + ,z − ) , the left side of Eq. (D8)increases from z − whereas the right side decreases from z + .3As a result, the gap between the two sides becomes smaller,and there cannot be any integer point in between.We now show that, for λ ( k ,k n ) < λ < λ ( z + ,z − ) , set P o is guaranteed to be nonempty. In this region of λ , wehave g ( k , k n ) > and g ( z + , z − ) < . Consider thepoint ( k a , k b ) = ( z + , k n ) . For g ( z + , k n ) = 0 , accordingto Eq. (D3), set P o is nonempty. The other two possibilities: g ( z + , k n ) > and g ( z + , k n ) < , can be treated separately.Suppose g ( z + , k n ) > , we can pick k i = z + , k i = z − and k i = k n . In this case, g ( k i , k i ) < and g ( k i , k i ) > hold by definition. It can then be shown that these two in-equalities imply g ( k i , k i ) < . In particular, note that g ( k i , k i ) − g ( k i , k i )= ( k i − k i ) (cid:0)(cid:0) λ − λ z (cid:1) k i − λz (cid:1) . (D18)As g ( k i , k i ) < , we have (cid:0)(cid:0) λ − λ z (cid:1) k i − λz (cid:1) k i < λzk i − z. (D19)Since k i = z − ≤ z and λz < , the right side is nega-tive and we have (cid:0)(cid:0) λ − λ z (cid:1) k i − λz (cid:1) < . This im-plies g ( k i , k i ) < g ( k i , k i ) < . For the other caseof g ( z + , k n ) < , we pick k i = k , k i = z + and k i = k n , so g ( k i , k i ) < and g ( k i , k i ) > hold bydefinition. From λz < and k i = z + ≥ z , we have (cid:0)(cid:0) λ − λ z (cid:1) k i − λz (cid:1) > . It can thus be concluded fromEq. (D18) that > g ( k i , k i ) > g ( k i , k i ) .The above proof procedure can be applied to the case ofpicking ( k i , k i , k i ) for λ ( k ,k n ) < λ < λ ( z + ,z − ) . Sup-pose there are four degrees k a > k b > z > k c > k d , with g ( k a , k d ) < and g ( k b , k c ) > . If the point ( k b , k d ) has g ( k b , k d ) > , we have that ( k i , k i , k i ) = ( k b , k c , k d ) de-fines a physical degree distribution from Eq. (33). Similarly,if g ( k a , k c ) < , we can choose ( k i , k i , k i ) = ( k a , k b , k c ) .The results of this appendix can be summarized as follows.Let λ = λ ( k ,k n ) and λ = λ ( z + ,z − ) . For λ < λ < λ , theset P o is empty and we can only find the local maxima amongall degree distributions that are supported on one or two nodaldegrees. For λ ≤ λ ≤ λ , set P o is nonempty, and it isnecessary to further analyze if there are other local maximasupported on one or two nodal degrees and if the degree dis-tributions in P o are maxima and global maxima. For λ > λ ,set P o again becomes empty. Appendix E: Solution of the HMF equation for degreedistribution supported on two nodal degrees
The equilibrium point x ∗ is given by the solution of λk i (1 − x ∗ i )Θ ∗ = x ∗ i , λk i (1 − x ∗ i )Θ ∗ = x ∗ i . (E1)Multiplying the two equations in Eq. (E1) by p i and p i , re-spectively, and summing them, we obtain ψ ( p ) = λz Θ ∗ − λz (Θ ∗ ) . (E2)To obtain ψ ( p ) , it suffices to find the value of Θ ∗ .Equation (E1) gives x ∗ i = λk i Θ ∗ λk i Θ ∗ , x ∗ i = λk i Θ ∗ λk i Θ ∗ . (E3)From the definition of Θ ∗ , we have z Θ ∗ = p i k i x ∗ i + p i k i x ∗ i . (E4)Substituting Eq. (E3) and the values p i and p i in Eq. (38)into Eq. (E4), we obtain the following quadratic equation for Θ ∗ : β (Θ ∗ ) + β Θ ∗ + β = 0 , (E5)with the coefficients, β = λ zk i k i ,β = λz ( k i + k i − λk i k i ) ,β = z − λzk i − λzk i + λk i k i . (E6)Noting that the second moment of the degree distribution is h k i = p i k i + p i k i = z ( k i + k i ) − k i k i , (E7)We have β = z − λ h k i < as λ is above the epidemicoutbreak threshold z/ h k i . Since β > , the only physicalsolution (with < Θ ∗ < ) of Eq. (E5) is Θ ∗ = − β + p β − β β β . (E8)Substituting Eq. (E8) into Eq. (E2), we obtain the value of ψ ( p ) as given by Eq. (39). [1] J. Aguirre, D. Papo, and J. M. Buld´u, “Successful strategies forcompeting networks,” Nat. Phys. , 230 (2013).[2] L. Pan, W. Wang, S. Cai, and T. Zhou, “Optimal interlayerstructure for promoting spreading of the susceptible-infected-susceptible model in two-layer networks,” Phys. Rev. E ,022316 (2019).[3] R. Pastor-Satorras, C. Castellano, P. Van Mieghem, andA. Vespignani, “Epidemic processes in complex networks,”Rev. Mod. Phys. , 925 (2015). [4] C. Castellano, S. Fortunato, and V. Loreto, “Statistical physicsof social dynamics,” Rev. Mod. Phys. , 591 (2009).[5] R. Pastor-Satorras and A. Vespignani, “Epidemic dynamics andendemic states in complex networks,” Phys. Rev. E , 066117(2001).[6] R. Cohen, S. Havlin, and D. ben-Avraham, “Efficient im-munization strategies for computer networks and populations,”Phys. Rev. Lett. , 247901 (2003).[7] R. Pastor-Satorras and A. Vespignani, “Immunization of com- plex networks,” Phys. Rev. E , 036104 (2002).[8] M. Kitsak, L. K. Gallos, S. Havlin, F. Liljeros, L. Muchnik,H. E. Stanley, and H. A. Makse, “Identification of influentialspreaders in complex networks,” Nat. Phys. , 888 (2010).[9] L. L¨u, D. Chen, X.-L. Ren, Q.-M. Zhang, Y.-C. Zhang, andT. Zhou, “Vital nodes identification in complex networks,”Phys. Rep. , 1 (2016).[10] T. W. Valente, “Network interventions,” Science , 49 (2012).[11] S. Goel, A. Anderson, J. Hofman, and D. J. Watts, “The struc-tural virality of online diffusion,” Management Sci. , 180(2016).[12] E. Ferrara, O. Varol, C. Davis, F. Menczer, and A. Flammini, “The rise of social bots,” Commun. ACM , 96 (2016).[13] S. Vosoughi, D. Roy, and S. Aral, “The spread of true and falsenews online,” Science , 1146 (2018).[14] C. Musco, C. Musco, and C. E. Tsourakakis, “Minimizing po-larization and disagreement in social networks,” in Proceedingsof the 2018 World Wide Web Conference (2018) pp. 369–378.[15] L. Wang and G.-Z. Dai, “Global stability of virus spreading incomplex heterogeneous networks,” SIAM J. Appl. Math. ,1495 (2008).[16] G. St-Onge, J.-G. Young, E. Laurence, C. Murphy, and L. J.Dub´e, “Phase transition of the susceptible-infected-susceptibledynamics on time-varying configuration model networks,”Phys. Rev. E97