Impacts of clustering on interacting epidemics

Journal of Theoretical Biology 304 (2012) 121–130

Contents lists available at SciVerse ScienceDirect

Journal of Theoretical Biology journal homepage: www.elsevier.com/locate/yjtbi

Impacts of clustering on interacting epidemics Bing Wang a,n, Lang Cao a, Hideyuki Suzuki a,b, Kazuyuki Aihara a a b

Institute of Industrial Science, The University of Tokyo, 4-6-1 Komaba, Meguro-ku, Tokyo 153-8505, Japan PRESTO, Japan Science and Technology Agency (JST), 4-1-8 Honcho Kawaguchi, Saitama 332-0012, Japan

a r t i c l e i n f o

abstract

Article history: Received 17 May 2011 Received in revised form 8 March 2012 Accepted 12 March 2012 Available online 28 March 2012

Since community structures in real networks play a major role for the epidemic spread, we therefore explore two interacting diseases spreading in networks with community structures. As a network model with community structures, we propose a random clique network model composed of different orders of cliques. We further assume that each disease spreads only through one type of cliques; this assumption corresponds to the issue that two diseases spread inside communities and outside them. Considering the relationship between the susceptible–infected–recovered (SIR) model and the bond percolation theory, we apply this theory to clique random networks under the assumption that the occupation probability is clique-type dependent, which is consistent with the observation that infection rates inside a community and outside it are different, and obtain a number of statistical properties for this model. Two interacting diseases that compete the same hosts are also investigated, which leads to a natural generalization of analyzing an arbitrary number of infectious diseases. For two-disease dynamics, the clustering effect is hypersensitive to the cohesiveness and concentration of cliques; this illustrates the impacts of clustering and the composition of subgraphs in networks on epidemic behavior. The analysis of coexistence/bistability regions provides significant insight into the relationship between the network structure and the potential epidemic prevalence. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Epidemiology Clustering coefficient Bond percolation Complex networks

1. Introduction The development of network methodology has led to the study of infectious diseases spread on contact networks. Unlike homogenous mixing, the incorporation of individual contacts induces new focuses on the problem of how network structures, such as clustering, influence the epidemic properties. Recently, network epidemiology has also been applied to cultural transmission processes in social systems, such as the spread of messages, rumors, or opinions, which have mechanisms qualitatively similar to the spread of diseases in some respects (Moreno et al., 2004; Zhou et al., 2007; Trpevski et al., 2010; Wang et al., 2012; Souza and Goncalves). Real networks are often characterized by communities in which vertices are more densely connected to each other than to the rest of the network (Girvan and Newman, 2002; Newman, 2004, 2006; Newman and Girvan, 2004; Reichardt and Bornholdt, 2006; Fortunato, 2010). Usually these communities are not completely separated but shared some vertices or edges (Palla et al., 2005). Vertices bridging communities play a major role in spreading dynamics (Salathe´ and Jones, 2010), because they can propagate different diseases/opinions through links inside communities or links that bridge communities; e.g., in online social

n

Corresponding author. Tel.: þ81 3 5452 6697. E-mail address: [email protected] (B. Wang).

0022-5193/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.jtbi.2012.03.022

network sites, such as forums, one user can send public messages in a general forum community, or one can send private messages to a closer friend. The spread of competing opinions/diseases has been broadly observed in social and biological networks (Wang et al., 2012; Souza and Goncalves; Huang and Rohani, 2006, 2005; Newman, 2005; Funk and Jansen, 2010; Danon et al., 2011; Karrer and Newman, 2011), where opinions can often be exclusive to each other, e.g., agreement or disagreement; and leprosy and tuberculosis are, for example, diseases such that if people are exposed to one of them, then they may get immune to the other (Karlen, 1996); strains of an infectious disease, e.g., influenza, have often ˜ o et al., 2005). cross-immunity (Nun Motivated by the above observations, we investigate two interacting diseases spreading on networks with community structures. As a network model with community structures, we propose a random clique network model described by a joint clique degree distribution, which is composed of different orders of cliques. It is natural to regard unions of cliques as communities as in the definition of the k-clique community proposed by Palla et al. (2005); each k-clique community is a union of k-cliques, which connect each other by rolling a k-clique through a series of adjacent k-cliques that share k1 vertices (Palla et al., 2005). For a network composed of m-cliques and n-cliques (m 4 n), by rolling an m-clique, m-cliques are detected as m-clique communities while n-cliques form links that connect different communities.

122

B. Wang et al. / Journal of Theoretical Biology 304 (2012) 121–130

These m-cliques overlap at vertices that connect them (Palla et al., 2005). This network model is the simplest one to characterize the community feature and helpful for theoretical analysis. For theoretical analysis based on the bond percolation, we further assume that each disease spreads only through one type of cliques for simplicity; this assumption represents the issue that two diseases spread inside communities and outside them through vertices that are overlapped by different communities. Even under this assumption, our model captures two important aspects of our original motivation: (1) competition of two diseases and (2) effect of community structures. What we are concerned with are the final outcome of the epidemics and the epidemic threshold rather than the precise time evolution, thus it is possible to map the problem of epidemic spreading onto a bond percolation problem, in such a way that the disease is transmitted along the occupied edges with the occupation probability being the average transmissibility; the ultimate size of the outbreak would be the size of a cluster being reached from an initial vertex by traversing the occupied edges (Newman, 2002, 2003, 2009; Kiss and Green, 2008; Gleeson, 2009; Miller, 2009a; Allard et al., 2009). Note that although there might be some imperfections for the application of bond percolation theory to the study of epidemics in the form of a fixed infectious period, which was then discussed and improved (Kenah and Robins, 2007; Miller, 2007, 2009b; Trapman, 2007), it still provides a simple and efficient way to understand the epidemic process on complex networks theoretically. In this work, we first generalize the clique random network model to include arbitrary order of cliques (Newman, 2009), which allows us to easily tune the degree distribution and the clustering coefficient. Based on the probability generating functions (pgfs) for the degree distribution and the excess degree distribution, we mathematically analyze different statistical properties including the size of the giant component, the distribution of component sizes, and the phase transition where the giant component first appears. We then study the epidemic behavior on this network model with bond percolation theory under the assumption that the occupation probability is clique-type dependent, which may be related to the observation that the infection rates inside of a community and outside of it are usually different. Based on these results, we investigate the interacting epidemics on clique random networks. We explore a possible existence region for each disease and find that in clique random networks, clustering as well as degree variance always decreases the critical threshold of the corresponding disease. However, when the average connectivity is fixed, spreading becomes more difficult in higher-order clique networks than in lower-order clique networks. The degree variance also produces a stable coexistence region for two diseases. The possibility of the final state falling into a coexistence/bistability region is analyzed by extending the definition of betweenness centrality to clique centrality. With the help of a clique centrality matrix, the relationship between network structure and possible interacting epidemics is also discussed.

2. The model Similarly to the definition of a degree, which is simply the number of links that a vertex has to other vertices, each vertex in a network has a ‘clique degrees’ from order 2 up to order M, denoted as fc2 , . . . ,cM g, each element of which cj is the number of j-clique that the vertex belongs to (see Fig. 1). Then pðc2 , . . . ,cM Þ means the fraction of vertices which have the specified arrangement of clique degrees fc2 , . . . ,cM g. A clique is defined as a complete subgraph in networks and its order is the number of

Fig. 1. Schematic illustration of the clique degree sequence of a vertex. Here, we do not consider 1-clique, which is composed of only one vertex. Since the network size is assumed to be infinitely large, the number of overlapping edges among different cliques is small and could be neglected. In this example, the clique degree sequence of the black vertex is fc2 ,c3 ,c4 g ¼ f1; 3,1g.

vertices contained in the subgraph. We define a j-clique as a clique where maximal order is j. Later in this paper, we occasionally mention a 3-clique/4-clique as a triangle/complete square, respectively, and the order of the clique that a link belongs to is also deemed to be its link type. The pgf for producing the clique degree distribution of a network is a function of M  1 variables (Wilf, 1994) as follows: g 0 ðxÞ ¼

1 X

pðc2 , . . . ,cM Þ

c2 ,...,cM

M Y

xcl l ,

ð1Þ

l¼2

where x ¼ fx2 , . . . ,xM g. Noting that each vertex in a j-clique corresponds to j1 links connecting to other members of the clique, the total degree for each vertex is k ¼ c2 þ 2c3 þ    þ ðM1ÞcM . The probability of finding a vertex with degree k in the network, can be written as pk ¼

1 X c2 ,...,cM

pðc2 , . . . ,cM Þdk, PM

j ¼ 2

ðj1Þcj

,

ð2Þ

where d:;: is the Kronecker delta. The pgf for the total degree distribution pk is 1 1 PM X X Pk ðzÞ ¼ pk zk ¼ pðc2 , . . . ,cM Þz l ¼ 2 ðl1Þcl ¼ g 0 ðz,z2 , . . . ,zM1 Þ: k

c2 ,...,cM

ð3Þ The clustering coefficient C measures the level of clustering in a network, and is defined as the average probability that two vertices are connected, given that they share a common neighbor in a network (Watts and Strogatz, 1998). Mathematically, this can be calculated as three times the ratio of the number of triangles Nn in the network to the number of connected triples of vertices N3 as follows: C¼

3Nn : N3

ð4Þ

Here, N n and N3 are      1 M M X X X l 1 l 1 @g 0  c pðc2 , . . . ,cM Þ ¼ n 3 3 3Nn ¼ n  3 l l 3 l @xl  c2 ,...,cM l ¼ 3 l¼3

,

x¼1

ð5Þ  1   X k n @2 P k ðzÞ pk ¼ N3 ¼ n  2 2 @z2  k

, z¼1

ð6Þ

B. Wang et al. / Journal of Theoretical Biology 304 (2012) 121–130

where n is the total number of vertices, nð1=lÞ@g 0 =@xl 9x ¼ 1 represents the total number of l-cliques (divided by l because there are l nodes in each l-clique), and ð3l Þ means the number of triangles in each l-clique. Substituting Eqs. (5) and (6) into Eq. (4), we obtain    l 1 @g 0  PM  l¼33 3 l @xl x ¼ 1 : ð7Þ C¼  1 @2 Pk ðzÞ  2 2 @z z¼1 Another important quantity we will use is the excess clique degree distribution, which is formally defined as the degree distribution of a vertex reached by following a randomly chosen edge. In the present study, there are totally M 1 types of cliques that could be followed,1 e.g., single links, triangles, or complete squares. We therefore need to define the excess clique degree distributions for each type of cliques. Assume qj ðc2 , . . . ,cM Þ, j Z2 is the distribution of the number of single links, triangles, etc., attached to a vertex reached by following a random j-clique; excluding the traversed clique, we have qj ðc2 , . . . ,cM Þ ¼

ðcj þ 1Þpðc2 þ d2,j , . . . ,cM þ dM,j Þ /cj S

ðj Z2Þ,

ð8Þ

where /cj S is the average number of j-cliques. Then the pgf to produce this distribution is defined as g j ðxÞ ¼

1 X

qj ðc2 , . . . ,cM Þ

c2 ,...,cM

M Y

xcl l ¼

l¼2

1 @g 0 /cj S @xj

ðj Z 2Þ:

ð9Þ

M1

g 0 ðf 2 , . . . ,f M

1 X

Þ¼

123 M Y

pðc2 , . . . ,cM Þ

c2 ,...,cM

ðf j Þðj1Þcj :

ð13Þ

j¼2

The expected size S of the giant component as a fraction of the entire network is M1

S ¼ 1g 0 ðf 2 , . . . ,f M

Þ,

ð14Þ

where ff 2 , . . . ,f M g are the solutions of Eqs. (10)–(12), which can be obtained by fixed-point iteration with appropriate initial values. 2.2. The average size of components To determine the expected size of connected components in the graph, we assume that hj(z) is the pgf for the distribution of the sizes of components that are reached by following a random j-clique (j Z2). Since hj(z) excludes the giant component, before the appearance of the giant component, the small component can be taken locally as a tree-like structure at the level of a single link, a triangle, or a higher-order clique, and therefore, along a random clique a vertex is reached, by which the cluster extends to any number of other clusters generated with the same distribution. For each type of cliques being traversed, the pgf hj(z) satisfies the following condition: M1

hj ðzÞ ¼ zg j ðh2 ðzÞ, . . . ,hM

ðzÞÞ,

ð15Þ

and the probability that a randomly chosen vertex belongs to a component of a given size is generated by the function h0 ðzÞ, M1

h0 ðzÞ ¼ zg 0 ðh2 ðzÞ, . . . ,hM

ðzÞÞ:

ð16Þ

Now we get a total of M 1 pgfs for the excess clique degree distributions. With the pgf for the clique degree distribution defined in Eq. (1), in the next section, we demonstrate the size of the giant component.

The average size of small components can be obtained by calculating the first derivative of h0 ðzÞ in Eq. (16), evaluated at z¼1, i.e., h00 ð1Þ,

2.1. The size of the giant component

h00 ð1Þ ¼ 1 þ

Consider the probability that a vertex is not connected to the giant component through a random j-clique link, denoted by fj; the probability that a j-clique does not belong to the giant

where g 0ð:;...;:;;...;:Þ represents the differentiation with respect to its arguments. h0j ð1Þ can be obtained by differentiating Eq. (15) with respect to z, evaluated at z ¼1 as follows:

j1

component is f j

. Note that a vertex separating from the giant

the clique degree distribution qj ðc2 , . . . ,cM Þ as the probability of reaching a vertex with clique degrees fc2 , . . . ,cM g along a random j-clique. We therefore have the following equations: 1 X

q2 ðc2 , . . . ,cM Þ

c2 ,...,cM

f3 ¼

1 X c2 ,...,cM

M Y

ðj1Þcj

¼ g 2 ðf 2 ,f 3 , . . . ,f M

ðj1Þcj

¼ g 3 ðf 2 ,f 3 , . . . ,f M

fj

2

M1

2

M1

Þ,

ð10Þ

Þ,

ð11Þ

j¼2

q3 ðc2 , . . . ,cM Þ

M Y

fj

j¼2

^

fM ¼

1 X c2 ,...,cM

qM ðc2 , . . . ,cM Þ

M Y

ðj1Þcj

fj

ðd ,..., dm,j ,..., dM,j Þ 0 hj ð1Þ,

ðj1Þg 0 2,j

ð17Þ

j¼2

component implies that all its connected neighbors do not belong to the giant component, and thus the probability that a vertex is Q ðj1Þcj not connected to the giant component is M . Consider j ¼ 2 fj

f2 ¼

M X

2

M1

¼ g M ðf 2 ,f 3 , . . . ,f M

Þ:

ð12Þ

h0j ð1Þ ¼ 1 þ

M X l¼2

ðl1Þ

Hjl 0 h ð1Þ, /cj S l

ð18Þ

where Hjl are the elements of the Hessian matrix H of the second derivatives of g0, evaluated at z ¼ 1. Then in a matrix form, this can be written as h ¼ 1þ a1 Hb  h, where h ¼ fh02 ð1Þ, . . . ,h0M ð1Þg> , 1 ¼ f1, . . . ,1g> , and a ¼ diagð/c2 S, . . . ,/cM SÞ, b ¼ diagð1, . . . ,M1Þ, 0 ð2;0,...,0Þ 1  g0 g ð1;1,0,...,0Þ    g ð1;0,...,1Þ  0 0 B ð1;1,0,...,0Þ C ð0;2,0,...,0Þ ð0;1,...,1Þ Bg C g    g B C 0 0 H¼B 0 : C B C ^ ^ ^ ^ @ A  g ð0;1,0,...,1Þ    g ð0;0,0,...,2Þ g ð1;0,...,1Þ  0 0 0 z¼1

The average size of a component can be calculated by solving Eq. (18) and substituting it into Eq. (17). The average size of a component is a finite value and it will diverge at the point satisfying the condition detðIa1 HbÞ ¼ 0, which corresponds to the appearance of a giant component.

j¼2

Then the probability that a randomly chosen vertex is not in the giant component is 1 In the present case of multiple clique types, the excess clique degree distribution is defined in a similar way as the clique degree distribution of a vertex reached by following a random clique, which should be proportional to the number of this type of cliques and normalized by the average number of them.

3. Bond percolation on clique random networks The problem of epidemic spread described by an SIR model can be mapped onto a bond percolation problem in such a way that the disease is transmitted along the occupied edges with probability T, which is the average transmissibility and the size of the outbreak is equivalent to the component reached by traversing

124

B. Wang et al. / Journal of Theoretical Biology 304 (2012) 121–130

the occupied edges from the initial infected vertex. Since the infection rates inside of a subgroup and outside of a subgroup are usually different, we therefore assume that the occupation probability is clique-type dependent and described in a vector form T ¼ ðT 2 , . . . ,T M Þ> . A randomly chosen vertex may connect to a giant cluster through neighboring vertices located at the ends of the occupied links, which are also in a giant cluster. For example, in a 3-clique (triangle), one corner vertex (the focused vertex in red) can be connected to a giant cluster via the neighbors of an occupied edge with probability T3 (Fig. 2(a)). It is also possible that the vertex at the end could be connected to a giant cluster via an occupied single link (a clique of one order lower than the original clique), and then the focused vertex will also be in the giant cluster (Fig. 2(b)). Similarly to the discussion in Section 2 of the giant component, we assume that the probability that a random vertex is not connected to a giant cluster via a single link is f2. If the vertex is not connected to the giant cluster, then either the edge is not occupied with probability 1T 2 or although the edge is occupied with probability T2, the vertex at the end of the link is not connected to the giant cluster by any other cliques via single links, triangles, etc., with the corresponding probability 2

M1

2

M1

g 2 ðf 2 ,f 3 , . . . ,f M Þ; for simplicity, we denote f ¼ ff 2 ,f 3 , . . . ,f M Therefore, we have f 2 ¼ 1T 2 ð1g 2 ðf ÞÞ:

g.

the end of the occupied edges are connected to a giant cluster with probability 1g jn ðf Þ, and with probability g jn ðf Þ, they are not connected to a giant cluster but they will be connected through their neighbors in cliques of one order lower with probability n1j

1f nj

ðg n ðf Þ; T n Þ.

For example, we can naturally derive the probability that a random vertex is not connected to a giant cluster via a 4-cliquetype link, f4, satisfying the condition, 3

f 4 ¼ 13T 4 ð1T 4 Þ2 ½1g 4 ðf Þ þ g 4 ðf Þyðf 4 ; T 4 Þ 3T 24 ð1T 4 Þf1g 24 ðf Þ þg 24 ðf Þ½T 4 ð1g 4 ðf ÞÞgT 34 ð1g 34 ðf ÞÞ,

ð22Þ

where

yðf 4 ; T 4 Þ ¼ 2T 4 ð1T 4 Þ½1g 4 ðf Þ þ g 4 ðf ÞT 4 ð1g 4 ðf ÞÞ þ T 24 ð1g 24 ðf ÞÞ is the probability that the vertex is connected to the giant cluster through neighboring vertices in 3-cliques. By solving Eq. (21), we can calculate the size of a giant cluster S as 2

M1

S ¼ 1g 0 ðf 2 ,f 3 , . . . ,f M

Þ:

ð23Þ

ð19Þ

Now we consider the probability that a random vertex is not connected to a giant cluster via a 3-clique-type link, denoted by f3. Then the probability that a vertex is connected to a giant cluster 2 by any of the corner vertices in a 3-clique is 1f 3 . Two cases must happen, and at least one link is occupied. We first consider the first case that there is only one link being occupied (see Fig. 2(a) and (b)). If the end point of the occupied edge is connected to a giant cluster (Fig. 2(a)), the vertex will be in the giant cluster with probability 2T 3 ð1T 3 Þð1g 3 ðf ÞÞ. The end point may be connected to a giant cluster only through its neighbors in cliques of one order lower, i.e., single links (Fig. 2(b)). Then the focused vertex will also be in the giant cluster with probability 2T 3 ð1T 3 Þg 3 ðf ÞT 3 ð1g 3 ðf ÞÞ. Another case is that two edges are both occupied, then the focused vertex will be in the giant cluster with probability T 23 ð1g 23 ðf ÞÞ (see Fig. 2(c) and (d)). Then we have 2

f 3 ¼ 12T 3 ð1T 3 Þ½1g 3 ðf Þ þ g 3 ðf ÞT 3 ð1g 3 ðf ÞÞT 23 ð1g 23 ðf ÞÞ:

ð20Þ

Based on the above results, we recursively derive the probability fn that a random vertex is not connected to a giant cluster via an n-clique-type link (n Z 2) by,  n1  X n1 j n1 T n ð1T n Þn1j ½ð1g jn ðf ÞÞ f n ðg n ðf Þ; T n Þ ¼ 1 j j¼1 n1j

þ g jn ðf Þð1f nj

ðg n ðf Þ; T n ÞÞ,

ð21Þ

0

with f 1 ¼ 1. For a vertex in an n-clique, j links among n  1 links may be occupied and the remainder is unoccupied. The vertices at

T3

1-T3

T3

1-T3

4. Applications to interacting epidemics on clique networks 4.1. Two interacting epidemics In order to investigate the issue that two diseases which compete the same hosts in the population spread inside of subgroups and outside of them, we explore this problem on clique random networks in a way that two diseases spread along different types of cliques. To study the effects of cohesiveness resulted by different orders of cliques, networks are composed of cliques up to the fourth order (single links and triangles/complete-squares, respectively). Based on the analysis of two interacting diseases, we generalize our model to include an arbitrary number of interacting diseases. The critical threshold is often described by the basic reproduction number R0, which is the expected number of secondary cases produced in a completely susceptible population, by a typical infective individual. If R0 o 1, then the disease cannot invade the population, but if R0 41, then invasion is always possible. In the following we focus on how to calculate the critical threshold for each disease. Given a network composed of two orders of cliques, for example, single links and 4-cliques, with an arbitrary clique degree distribution pds, the pgf for the clique degree distribution is X pds xd ys : ð24Þ g 0 ðx,yÞ ¼ ds

T3

T3

T3

T3

T3 Fig. 2. Schematic illustration of occupied modes in a 3-clique (triangle). Red represents the focused vertex, gray represents a vertex that is connected to a giant cluster by its neighbors in lower-order cliques, black represents a vertex in the giant cluster and white represents the vertex that is not connected to a giant cluster. Dark blue links represent the occupied links. In (a) and (b), only one link is occupied and the other link is unoccupied. (a) The vertex at the end of the occupied link is connected to a giant cluster. (b) The vertex at the end of the occupied link is connected to a giant cluster by a neighbor in a single-edge. In (c) and (d), both links are occupied and at least one or two neighbors are connected to a giant cluster. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

B. Wang et al. / Journal of Theoretical Biology 304 (2012) 121–130

The excess clique degree distributions of a vertex reached by following a random d-clique and s-clique are generated by P dp xd1 ys , ð25Þ g d ðx,yÞ ¼ dsP ds ds dpds P g s ðx,yÞ ¼

spds x dsP

d s1

y

ds spds

,

ð26Þ

respectively. Assume that the transmissibilities for the two diseases (we call them disease-I and disease-II) spreading along single links and 4-cliques are T2 and T4 and a host becomes immune to the other disease if he (she) is infected by one disease. Like the previous analysis, we consider the probability that a random vertex is not infected by a particular neighbor along a single link, f2, and the probability that a random vertex is not 3 infected by a 4-clique-type link, f4, then f 4 is the probability that the vertex is not infected by any corner vertices in a 4-clique. f2 and f4 can be solved in a similar way as shown in Eqs. (19) and (22). For two diseases, there should be enough susceptible individuals in the network, therefore, we have to consider a random vertex’s neighbors in the remaining networks. It should be noted that although we call it ‘remaining networks’, diseases spread simultaneously. The probability b2 , that its neighbors are not infected by disease-II either is given as follows: P dp ðf Þ3s P ds 4 : b2 ¼ ds ð27Þ ds dpds 3 f4

in the pgf of the excess clique degree function Replacing y with g d ðx,yÞ and applying an appropriate normalization factor, we obtain the pgf for the excess clique degree distribution of the vertices being not infected by disease-II, P 3s d1 3 ds dpds ðf 4 Þ x Gd ðx,f 4 Þ ¼ P : ð28Þ 3s ds dpds ðf 4 Þ If the edges in the remaining network are occupied with 3 probability T2, we need the pgf Gd ðx,f 4 ; T 2 Þ for the degree distribution of a vertex in the remaining network; this can be obtained by replacing x with 1T 2 þ T 2 x. The neighboring vertices should not be infected either and survive with probability b2 . The epidemic threshold can therefore be evaluated by setting the 3 condition b2 G0d ð1,f 4 ; T 2 Þ ¼ 1 (Newman, 2005), which can be simply understood as along a random single link a vertex who survives from disease-II is reached, and the disease passes through it and spreads out to infect at least one neighboring vertex who also survives from disease-II. Then the critical threshold of T2c for disease-I satisfies the following equation: P dðd1Þpds ðf 4 Þ3s b2 T 2 dsP ¼ 1: ð29Þ 3s ds dpds ðf 4 Þ Similarly, the probability that a random vertex is not infected by disease-I is f2, and the probability that its neighbors are not P P infected either is b4 ¼ ds spds ðf 2 Þd = ds spds , so the probability that 3

a 4-clique is not infected by disease-I is b4 . The pgf for the excess clique degree distribution of the vertex reached by following a random 4-clique is P d s1 ds spds ðf 2 Þ y : ð30Þ Gs ðf 2 ,yÞ ¼ P d ds spds ðf 2 Þ Let fs ðy; T 4 Þ be the function of the probability y that a vertex is not infected by its neighbors in a 4-clique, then

fs ðy; T 4 Þ ¼ 13T 4 ð1T 4 Þ2 f1yþ yyðy; T 4 Þg 3T 24 ð1T 4 Þ½1y2 þ y2 ðT 4 ð1yÞÞT 34 ð1y3 Þ,

yðy; T 4 Þ ¼ 2T 4 ð1T 4 Þ½1y þ yT 4 ð1yÞ þT 24 ð1y2 Þ

ð31Þ

where is, as used in Eq. (22), the probability that a vertex is connected to the giant

125

cluster through its neighboring vertices in 3-cliques. By setting 3

the condition b4 G0s ðf 2 ,1; T 4 Þ ¼ 1, we then obtain the critical transmissibility T4c satisfying the following condition: P sðs1Þpds ðf 2 Þd b34 f0s ð1; T 4 Þ dsP ¼ 1: ð32Þ d ds spds ðf 2 Þ By solving Eqs. (29) and (32), we obtain the critical threshold for each disease. 4.2. Generalization to an arbitrary number of interacting epidemics The above analysis can be generalized to the interaction of an arbitrary number of diseases on clique random networks, each of which spreads along a kind of clique up to order M. Assume that the joint clique degree distribution follows pðc2 , . . . ,cM Þ, the transmissibility for an n-clique is Tn, and each host becomes immune to other kinds of diseases once one is infected by a specific disease. To calculate the critical threshold for each j1 disease, we first calculate the probability f j that a random vertex is not infected by a j-type disease. Based on the previous results, the pgf for the effective excess clique degree distribution for an n-clique is defined as QM P ðl1Þcl cn xn l ¼ 2,l a n ðf l Þ c ,...,c cn pðc2 , . . . ,cM Þ M1 g n ðf 2 , . . . ,xn ,f M Þ ¼ P2 M : QM ðl1Þcl c pðc , . . . ,c Þ ðf Þ n M 2 l l ¼ 2,l a n c2 ,...,cM ð33Þ To calculate the critical transmissibility we have to derive the probability that a neighbor vertex is not infected by other diseases either, bn , as follows: Q P ðl1Þcl cn pðc2 , . . . ,cM Þ M l ¼ 2,l a n ðf l Þ : ð34Þ bn ¼ c2 ,...,cM P c2 ,...,cM cn pðc2 , . . . ,cM Þ The critical transmissibility for disease spreading along an ncliques, Tnc is then given by QM P ðl1Þcl l ¼ 2,l a n ðf l Þ c2 ,...,cM cn ðc n 1Þpðc2 , . . . ,c M Þ 0 bn1 f ð1; T Þ ¼ 1, QM P n n ðl1Þcl l ¼ 2,l a n ðf l Þ c2 ,...,cM cn pðc2 , . . . ,cM Þ ð35Þ where similarly to Eq. (21), fn ðxn ; T n Þ is defined as  n1  X n1 l T n ð1T n Þn1l fn ðxn ; T n Þ ¼ 1 l l¼1 ½1xln þ xln ð1fnl ðxn ; T n ÞÞ,

ð36Þ

with f1 ¼ 1. 4.3. Numerical results for two interacting epidemics Clustering effect is naturally obtained by the introduction of different orders of cliques (cohesiveness), described by a joint clique degree distribution. This leads us to seek a way generating two kinds of clique networks incorporating clustering and degree heterogeneity, which allows the analysis on how clique cohesiveness, the clustering degree, and heterogeneity, take effect on the epidemic behavior of two interacting diseases. Let us start with two clique networks composed of single links and 3-cliques/4-cliques. There are two reasons for choosing theses networks. Firstly, two diseases spread inside of a community and outside of it can be equivalently investigated on clique networks where one disease spreads along higher order cliques like 3-cliques or 4-cliques, while the other disease spreads along lower order cliques like single links; secondly, 3-cliques and 4-cliques show different cohesiveness. Following Volz (2011), we then introduce a degree variance in networks with negative binomial (NB) distributions, which allows us to keep the mean

126

B. Wang et al. / Journal of Theoretical Biology 304 (2012) 121–130

1 2

3

0.6

1

0.4

2

0.6

1

0.4

0.2 0

C=0.05,ν=21 3

0.8 T3

0.8 T3

1

C=0.05,ν=18.5

0.2 0

0.2

0.6

0.4

0.8

0

1

0

0.2

1

T3

T3

0.2 0.6

0.4

3

0.2

1 0.2

2

0.6 0.4

3

0

1

C=0.10,ν=21

0.8

0.4

0

1

2

0.6

0.8

T2

C=0.10,ν=18.5

0.8

0.6

0.4

T2

0.8

1

T2

0

1 0

0.2

0.6

0.4

0.8

1

T2

Fig. 3. Comparison of possible existence regions for two interacting epidemics on networks composed of single links and 3-cliques with different C and n. T2 and T3 represent the occupation probabilities along single links and triangles. Each curve represents the critical threshold of one disease as a function of the occupation probability of the other disease. In the lower-left region surrounded by the horizontal and vertical lines, neither disease can survive. Regions denoted by ‘1’ and ‘2’ are the existence areas for disease-I and disease-II; region ‘3’ is an uncertain region, where each disease may exist or the two diseases coexist. The average connectivity is /kS ¼ 6.

value fixed while changing the variance. The negative binomial distribution NB(p, r) with parameters p and r is given by   kþ r1 r p ð1pÞk , pk ¼ ð37Þ r and the pgf for NB(p, r) is  r p g NB ðxÞ ¼ : 1ð1pÞx

ð38Þ

We can modify the distribution with a tunable fraction pt of edges constituting a triangle while keeping the mean of the distribution constant. The probability that a pair of edges appear as a triple instead of a triangle is 1pt . So for a random number k, 2-tuples are generated by Eq. (38), and the numbers of single links and triangles are generated by ðð1pt Þx2 þpt yÞk , where x is the variable for single edges and y is the variable for triangles. The pgf for the clique degree distribution of a network composed of single links and 3-cliques is g 3 ðx,yÞ ¼ g NB ðð1pt Þx2 þ pt yÞ. Similarly, the pgf for the clique degree distribution for a network composed of single links and 4-cliques is g 4 ðx,zÞ ¼ g NB ðð1ps Þx3 þ ps zÞ, where z is the variable for 4-cliques and ps is the probability that a triple of edges appear in a complete square. With given parameters p, r, and pt/ps, the properties of the network topology such as the average connectivity /kS, the degree variance n, and the clustering coefficient C can be easily derived by the pgf method. The parameters in this study are chosen such that the average connectivity is large enough to ensure that the network is connected.2 2 The mean of the NBðp,rÞ distribution is l ¼ rð1pÞ=p and the variance is v ¼ rð1pÞ=p2 , so for the network composed of single links and 3-cliques, the average connectivity /kS ¼ 2l, which can be calculated by the univariate pgf G3 ðxÞ ¼ g 3 ðx,x2 Þ, with /kS ¼ G03 ð1Þ, and the degree variance is n ¼ G3 00 ð1ÞG03 ð1Þ2 . Similarly, the average connectivity and the degree variance for the network composed of single links and 4-cliques can be calculated with the function

Eqs. (29) and (32) give the conditions that the critical transmissibility for each disease spreading on the network composed of 2-cliques and 4-cliques should satisfy, which are functions of the transmission rate of the other disease. In a similar way, the critical transmissibilities for two diseases spreading on networks composed of 2-cliques and 3-cliques can be obtained. In the following, the results are obtained by solving these equations. Fig. 3(a) and (c) compare the impacts of clustering on the existence regions of the two diseases. The curve indicates the critical transmissibility for each disease. The figure is divided into roughly four parts. The regions denoted by ‘1’ and ‘2’ represent the existence regions for the diseases which spread along single links (disease-I) and 3-cliques (disease-II), respectively. Note that the curves are horizontal and vertical for small T2 and T3. These constant lines correspond to the basic critical thresholds for diseases-II/I, which can be understood as circumstances under which the transmission rate of one disease is not large enough to influence the other and hence, it takes a constant value irrespective of the other one. With an increase in clustering (see Fig. 3(c)), more 3-cliques and fewer single edges appear in the network, as expected, inducing a lower basic critical threshold (horizontal line) for disease-II and a higher critical threshold (vertical line) for disease-I. After reaching the critical transmission value for each disease, the spread of one disease is influenced by the other by competing the common hosts of population in networks. Therefore, the critical transmissibility is a monotonically increasing curve as a function of the other disease’s transmission rate. The upper-right region between the two critical curves, denoted by ‘3’, is the uncertain region where each disease may exist or the two

(footnote continued) G4 ðxÞ ¼ g 4 ðx,x3 Þ. The clustering coefficient C for two networks can be obtained by ðzÞ ðÞ 00 00 calculating g ðyÞ 3 ð1; 1Þ, G3 ð1Þ, g 4 ð1; 1Þ, and G4 ð1Þ, respectively, as in Eq. (7). g 3;4 ð1; 1Þ means the first derivative of its argument evaluated at (1, 1).

B. Wang et al. / Journal of Theoretical Biology 304 (2012) 121–130

corresponding disease in each network. Although the high density in higher-order cliques speeds up the spreading of disease inside, the lack of candidate infection paths may prevent epidemic incidence for the corresponding disease. Heterogeneity in the degree distribution is confirmed to decrease the critical transmissibility and to introduce a stable coexistence region for the two diseases. With regard to the uncertain region of coexistence/ bistability, the percolation theory fails to predict which state the epidemic behavior will fall into. The comparison of epidemics on different orders of clique networks may also provide some insight for preventing the spread of diseases in reality. The efficiency of vaccine strategies can be usually measured by the basic reproduction ratio R0. We have found that although 4-cliques are denser than 3-cliques, the lack of enough infection routes (the number of 4-cliques should be less than that of 3-cliques because of the constraint of clustering and the total number of links in networks) halts successive disease spreading at some vertices. The intuitive notion for implementing an efficient vaccine strategy is to vaccinate vertices bridging different cliques. Clique networks are closely related to community networks, since the latter are relaxed forms of cliques in the sense of the density of internal edges to the number of external edges. They have often been observed in social and biological networks. The effect of community structure on disease dynamics has been studied in Salathe´ and Jones (2010), and it is suggested that immunization interventions targeted at individuals bridging communities (which can be measured by betweenness centrality) are more effective than targeting highly connected individuals. In Section 4.4, we will make use of the definition of vertex betweenness centrality and extend it to clique centrality, which will help the analysis of the relationship between the network structure and the epidemic behavior.

C=0.05,ν=18.5

2

T4

T4

diseases may coexist, depending on the network structure and the position of the initial infected vertex. We will discuss region ‘3’ in detail later. With increasing the degree variance, we observe a decreasing critical transmissibility for disease-II (see Fig. 3(b) and (d)), which is consistent with the previous results that the critical transmissibility vanishes in scale-free networks when the variance diverges to infinity (Newman, 2002). To provide some insight into the interactions between clustering and the composition of substructures in networks, we perform a similar calculation on networks with fixed connectivity but with the clique-type being replaced with 4-cliques. As shown in Fig. 4(a) and (c), the critical threshold curves show similar characteristics as in Fig. 3(a) and (c), i.e., clustering decreases the epidemic threshold for the corresponding disease. However, we observe a lower critical threshold for the disease spreading along single links, and a higher critical threshold in 4-clique networks than that in 3-clique networks. This can be understood as compared to 3-clique networks, 4-clique networks are denser and it is easier to get higher clustering, which releases a larger number of single links in networks due to the constraint of the constant number of total links. As a result, the critical threshold of the corresponding disease spreading along single links is decreased. Although the high density structure in 4-cliques is expected to speed up transmission, because of their limited number, a disease is still difficult to spread. With a further increase in clustering (see Fig. 4(e)), disease-II spreads more easily, with a decreased critical threshold, and the two threshold curves intersect and form a stable coexistence region ‘4’, where the two diseases can coexist in the network. The above results show that clustering and the composition of cliques in networks can have a substantial effect on the epidemic behavior. As expected, clustering decreases the critical transmissibility of the

1 0.8 0.6 0.4 0.2 0

1 0

0.2

0.4

0.6

0.8

1

1 0.8 0.6 0.4 0.2 0

2

C=0.05,ν=21

3

0

1

0.2

0.4

2

3

C=0.10,ν=18.5 1

0

0.2

0.4

0.6

0.8

1

1 0.8 0.6 0.4 0.2 0

2

3

1 0

0.2

0

0.6 T4

T4

2

0.2

4 1

0

0.2

0.4

0.6 T2

1

0.4

0.6

0.8

1

T2

C=0.21,ν=18.5

0.4

0.8

C=0.10,ν=21

T2 0.6

0.6 T2

T4

T4

T2 1 0.8 0.6 0.4 0.2 0

127

0.8

1

C=0.21,ν=21

0.4

2

0.2 0

0

0.2

0.4

3

4

1 0.6

0.8

1

T2

Fig. 4. Comparison of possible existence regions for two interacting epidemics on networks composed of single links and 4-cliques for different C and n. T2 and T4 represent the occupation probabilities along single links and complete squares. Each curve represents the critical threshold for each disease obtained by solving Eqs. (29) and (32). In the lower-left region surrounded by the horizontal and vertical lines, neither disease can survive. Regions denoted by ‘1’ and ‘2’ are the existence areas for disease-I and disease-II; region ‘3’ is an uncertain region, where each disease may exist or the two diseases coexist. Region ‘4’ indicates the region where the two diseases coexist. The average connectivity is /kS ¼ 6.

128

B. Wang et al. / Journal of Theoretical Biology 304 (2012) 121–130

4.4. Coexistence/bistability phase Next we concentrate our attention on the ‘coexistence/bistability’ region of the phase diagram, in which the percolation theory predicts the birth of a giant component, namely all vertices connecting to the giant component formed by the infectious vertices of each disease with a non-zero probability. However, there are two scenarios with respect to which type of infectious vertices constitute the giant component: (a) two pathogens can both spread throughout the network corresponding to the coexisting phase in epidemic dynamics; and (b) only one pathogen can spread out and ultimately dominate the whole population corresponding to the bistable phase. The coexistence/bistability threshold of transmissibilities for two competing pathogens from the percolation theory provides insight as to the role of clique structures in networks; these structures heavily affect the prevalence of these diseases. It is natural to take into account the infection paths for each disease, which trace from only a single initially infectious vertex in the network to all other susceptible vertices in its connected component by traversing the same clique-type links. Note that each vertex of clique degrees fc2 , . . . ,cM g lies on infecting paths for different M1 diseases. Also note the intuitive fact that non-overlapping infectious paths avail to the coexistence of two competing pathogens, and conversely, that mixed and crossed paths avail to the bistability. Therefore, the clique degree distribution, i.e., the clique structure, determines the extent to which infection paths cross, and hence into which phase region, i.e., coexistence or bistability, and which final state the epidemic dynamics falls. To better understand the effects of clique structure on the prevalence of competing diseases, we take an illustrative example that compares two very specific network topologies: (a) a multitype-star network and (b) a multitype-leaf network. Here we extend typical star-like networks to multitype-star/leaf networks by adding a new nodal type of duple star/leaf, which is composed of two connected vertices as a whole sharing the same other neighboring vertices, as shown in Fig. 5. Leaf vertices in a pair can only reach other leaf vertices in a pair through star vertices. In the extreme case of star-like networks, a typical path starts from an initially infectious leaf, then goes through a star, and ends up by reaching another leaf. One can easily verify that a multitype-star network has non-overlapping infectious paths for both pathogens which guarantee the coexisting phase, whereas all infectious paths in a multitype-leaf network sharing the common star are completely crossed, which leads to the bistable phase; the final

epidemic state depends on which disease ‘wins’ to percolate to the star. Learning from the concept of betweenness centrality, which describes the extent to which a vertex lies on the shortest paths between other vertices (Thadakamalla et al., 2005), we introduce clique centrality based on infectious paths for each disease. For each vertex k, consider the (local) (iþ1)-clique centrality bki which counts the fraction of infectious paths for the (i þ1)-clique/ disease. More precisely, k

bi ¼

X sst ðkÞ i , sst ðkÞ

ð39Þ

sakat s,t A local world

where sst i ðkÞ is the number of shortest paths along (iþ1)-cliquetype links from vertex s to t going through vertex k, and sst ðkÞ is the total number of shortest paths from vertex s to t going through vertex k. Here the local world is typically assumed to be the neighboring information of the root (the focused vertex k), i.e., the local world in the definition of local betweenness centrality adopts the set of a root vertex, and its first and second neighbors to lighten the burden of calculation (Thadakamalla et al., 2005). For further simplicity and convenience in applying the generating function method, we adopt the set of the first/ nearest neighbors of vertex k as its local world, and the average (iþ 1)-clique centrality /bii S can be approximated by /bii S ¼ P

/ickiþ 1  ðickiþ 1 1ÞSk k j /jc j þ 1

 ðjckjþ 1 1ÞSk þ sj a i /ickiþ 1  jckjþ 1 Sk

¼

hii

Sij hij

, ð40Þ

where /  Sk denotes the average over all vertices of the network, hij ¼ /ickiþ 1  jckjþ 1 Sk di,j /ickiþ 1 Sk , and ckiþ 1 denotes the number of (iþ 1)-cliques passing through vertex k. Then the denominator indicates the total number of links in all types of cliques passing through vertex k. Note that hij is the element of the Hessian Þ, evaluated matrix H of the second derivatives of g 0 ðx2 ,x23 , . . . ,xM1 M at x ¼ 1, where g 0 ðxÞ is the pgf of the clique degree distribution pðc2 , . . . ,cM Þ. Returning to the example in Fig. 5, given a multitype-star network with l leaves, we have its pgf of the clique degree l l distribution g star 0 ðx2 ,x3 Þ ¼ ðl=ðl þ 3ÞÞx2 x3 þð1=ðl þ 3ÞÞx2 þð2=ðl þ 3ÞÞx3 . Similarly, given a multitype-leaf network with l1 single leaves and l2 duple leaves, we have g leaf 0 ðx2 ,x3 Þ ¼ ðl1 =ðl1 þ 2l2 þ 1ÞÞx2 þð2l2 = ðl1 þ 2l2 þ 1ÞÞx3 þ ð1=ðl1 þ2l2 þ 1ÞÞxl21 xl32 . Therefore the Hessian

Fig. 5. Two illustrative network topologies with star/leaf-like structure throughout which two competing pathogens both spread (‘coexistence’), or one of these pathogens spreads out and ultimately dominates the whole population (‘bistability’), respectively. (a) In the multitype-star network for the ‘coexistence’ case, each leaf vertex () connects with a single star ðmÞ and duple stars ð~Þ; (b) in the multitype-leaf network for the ‘bistability’ case, a star vertex ðmÞ connects with all single leaves () and duple leaves ð’Þ. The lines represent different types of links contained in 2-cliques (thick) and 3-cliques (thin).

B. Wang et al. / Journal of Theoretical Biology 304 (2012) 121–130

By direct calculation of Hstar=leaf , we obtain the clique centrality for each network, 8 lðl1Þ > > if i ¼ j ¼ 1, > > > C > > < il þ i2 lðil1Þ star /bij S ¼ if i ¼ j 4 1, > C > > > > ijl > > if i a j, : C

matrices of the multitype-star/leaf network are given by 0

lðl1Þ B lþ3 star B H9x ¼ 1 ¼ B @ 2l lþ3

1 2l C lþ3 C C, 2l þ 4lð2l1Þ A lþ3

0

1 2l1 l2 l1 þ 2l2 þ 1 C C C, 4l2 þ 2l2 ð2l2 1Þ A l1 þ 2l2 þ 1

l1 ðl1 1Þ B l þ2l þ 1 2 leaf B 1 H9x ¼ 1 ¼ B @ 2l1 l2 l1 þ2l2 þ 1

star

ð41Þ

star

and thus we have /b11 S ¼ ðl1Þ=ð9l þ 1Þ, /b22 S ¼ ð8l2Þ= ð9l þ 1Þ for the multitype-star network, and /b11 Sleaf ¼ l1 ðl1 1Þ= 2 ððl1 þ 2l2 Þ2 l1 þ 2l2 Þ, /b22 Sleaf ¼ ð4l2 þ2l2 Þ=ððl1 þ 2l2 Þ2 l1 þ 2l2 Þ for the multitype-leaf network. As discussed above, the diagonal element /bii S or formally hnorm,ii of the normalized Hessian matrix Hnorm (all matrix elements sum up to one) describes the fraction of (iþ1)-clique infectious paths, which appears as an index of spreading coverage for the (iþ1)th disease, i.e., the extent to which the pathogen can spread out appreciably to cover a giant component. On the other hand, note that its off-diagonal elements /bij S are the fractions of such hybrid paths reaching vertex k along an (iþ1)-clique infectious path and then going out of vertex k along a (j þ1)clique infectious path, averaged over all vertices k in the network. Therefore these off-diagonal elements reflect the clique structural feature that avails to the emergence of a giant component for a single disease (corresponding to the bistable phase). As an example, consider the extension of multitype-star/leaf networks to an arbitrary clique order M, where the multitype-star extended network contains different types of M  1 stars(up to the M-clique) and l leaves connecting with all the stars, and the multitype-leaf network includes li pieces of i-polyploid leaves that are linked to the star with (iþ1)-cliques, for i ¼ 1, . . . ,M1. The pgfs of the clique degree distributions for multi-star network and multi-leaf network are g star 0 ðx2 ,x3 ,

l1 2l2 ðM1ÞlM1 x2 þ x3 þ    þ xM B B B 1 l1 l2 x x    xlMM1 , B 2 3

where A ¼ l þ MðM1Þ=2 and B ¼ star

of Hessian matrices H are 8 lðl1Þ > > if > > > A > > < 2 il þ i lðil1Þ star hij ¼ if > A > > > > ijl > > if : A

leaf

and H

PM1 ,

i ¼ j ¼ 1, i ¼ j 4 1, i aj,

and

leaf

hij

8 l1 ðl1 1Þ > > > > B > > > < 2 i li þ ili ðili 1Þ ¼ > B > > > > ili jlj > > : B

if i ¼ j ¼ 1, if i ¼ j 41, if ia j:

8 l1 ðl1 1Þ > > > > D > > > < 2 i li þ ili ðili 1Þ leaf /bij S ¼ > D > > > > il jl > i j > : D

if i ¼ j ¼ 1, if i ¼ j 41, if i aj,

where C¼

M 1 X

2

k lðkl1Þ þ

k¼1

M 2 X ðM þ 1ÞðM2Þ lþ kðM þ kÞðMk1Þl, 2 k¼1

and D¼

M 1 X

klk ðklk 1Þ þ

k¼1

M 1 X

2

k lk þ

k¼2

M 2 X k¼1

klk

M 1 X

mlm :

m ¼ kþ1

As the network size of the multitype-star approaches infinity (l-1), the corresponding Hstar norm approaches a diagonal matrix, star star /bii S  i3 ,/bi a j S-0, which produces the coexistence of each disease in the final epidemic state. For the multitype-leaf netleaf work, /bij S takes a non-ignorable value. This non-diagonally dominant Hleaf norm of the multitype-leaf network implies the clique structure of networks on which the epidemic dynamics tends to fall into a bistable phase. Moreover, combining with the occupation probability of each cliquey type link, we define the centrality matrix for infectious paths as ð44Þ

>

ð42Þ

þ

and

B ¼ Hnorm JTT> ¼ ð/bij ST i T j Þij ,

l 1 2 M1 l xM , . . . ,xM Þ ¼ x2 x3    xM þ xl2 þ xl3 þ    þ A A A A

g leaf 0 ðx2 ,x3 , . . . ,xM Þ ¼

129

where ðÞ and J denote the matrix transpose and Hadamard product (element-wise product), respectively. Although rigorous analysis is still required, whether the final state falls into the coexistence/bistability phase probably depends on the diagonally dominant feature of the path centrality matrix B.

ð43Þ

k ¼ 1 klk þ 1. Then the elements hstar and hleaf (for 1 r i,j rM1) ij ij

5. Conclusions Since community structures have been found to make a major influence on spreading dynamics, we explore two interacting diseases spreading in networks with communities. Based on the definition of the k-clique community, we proposed a random clique network model with the basic feature of community structures, which is composed of different orders of cliques. For theoretical analysis based on the bond percolation, we assumed that two diseases spread along only one type of cliques, which also represents the issue that two diseases spread inside communities and outside them. We apply the bond percolation theory to the analysis of a disease spreading on clique random networks, under the assumption that the occupation probability is clique-type dependent. We have demonstrated the applicability of bond percolation theory to the issue of two interacting epidemic diseases on this network model and explored possible existence regions for each disease. As expected, clustering decreases the critical threshold for the corresponding disease as well as the degree variation. By

130

B. Wang et al. / Journal of Theoretical Biology 304 (2012) 121–130

comparing two kinds of networks composed of different cohesive cliques, we have found that although higher cohesiveness in higher-order clique networks speeds up the epidemic spreading, the limited number of infection routes results in a higher critical threshold. The uncertain region of coexistence/bistability was investigated by focusing on the network structure to extend the definition of the betweenness centrality to the clique centrality of infected paths. With specific examples of multitype-star/leaf networks, we have demonstrated that whether the final state falls into coexistence or bistability is closely related to the characteristics of the path centrality matrix. It should be noted that in a recently published paper (Volz et al., 2011) the authors proposed the ODEs method with an edgebased compartmental model to predict the final epidemics exactly. Although the result with the original bond percolation theory gives a biased epidemic size, both methods find that the final prevalence usually decreases with increasing clustering (or variance) (Volz et al., 2011). Therefore, bond percolation theory is useful for the study of interacting epidemics. The results we obtained would give some insights on understanding how the clustering, the degree variance, and the cohesiveness of cliques interact for two diseases spreading in networks with communities.

Acknowledgments We would like to thank the anonymous reviewers for letting us know recently published important references Karrer and Newman (2011) and Volz et al. (2011). This research is supported by the Japan Society for the Promotion of Science (JSPS) through a Grant-in-Aid for JSPS Fellows (21-09275) and partially supported by the Aihara Project, the FIRST program from JSPS, initiated by CSTP. References ¨ P.-A., Dube´, L.J., Pourbohloul, B., 2009. Heterogeneous bond Allard, A., Noel, percolation on multitype networks with an application to epidemic dynamics. Phys. Rev. E 79 (3), 036113. Danon, L., Ford, A.P., House, T., Jewell, C.P., Keeling, M.J., Roberts, G.O., Ross, J.V., Vernon, M.C., 2011. Networks and the epidemiology of infections. Interdiscip. Perspect. Infect. Dis., 284909. Fortunato, S., 2010. Community detection in graphs. Phys. Rep. 486, 75–174. Funk, S., Jansen, V.A.A., 2010. Interacting epidemics on overlay networks. Phys. Rev. E 81 (3), 036118. Girvan, M., Newman, M.E.J., 2002. Community structure in social and biological networks. Proc. Natl. Acad. Sci. USA 192, 7821–7826. Gleeson, J.P., 2009. Bond percolation on a class of clustered random networks. Phys. Rev. E 80 (3), 036107.

Huang, Y., Rohani, P., 2005. The dynamical implications of disease interference: correlations and coexistence. Theor. Popul. Biol. 68 (3), 205–215. Huang, Y., Rohani, P., 2006. Age-structured effects and disease interference in childhood infections. Proc. R. Soc. Lond. B Biol. Sci. 273, 1229–1237. Karlen, A., 1996. Man and Microbes, 1st ed. Simon and Schuster, New York. Karrer, B., Newman, M.E.J., 2011. Competing epidemics on complex networks. Phys. Rev. E 84, 036106. Kenah, E., Robins, J.M., 2007. Second look at the spread of epidemics on networks. Phys. Rev. E 76 (3), 036113. Kiss, I.Z., Green, D.M., 2008. Comment on ‘‘Properties of highly clustered networks’’. Phys. Rev. E 78 (4), 048101. Miller, J.C., 2007. Epidemic size and probability in populations with heterogeneous infectivity and susceptibility. Phys. Rev. E 76 (1), 010101. Miller, J.C., 2009a. Percolation and epidemics in random clustered networks. Phys. Rev. E 80 (2), 020901. Miller, J.C., 2009b. Spread of infectious disease through clustered populations. J. R. Soc. Interface 6, 1121. Moreno, Y., Nekovee, M., Pacheco, A.F., 2004. Dynamics of rumor spreading in complex networks. Phys. Rev. E 69, 066130. Newman, M.E.J., 2002. Spread of epidemic disease on networks. Phys. Rev. E 66 (1), 016128. Newman, M.E.J., 2003. Properties of highly clustered networks. Phys. Rev. E 68 (2), 026121. Newman, M.E.J., 2004. Detecting community structure in networks. Eur. Phys. J. B 38, 321–330. Newman, M.E.J., 2005. Threshold effects for two pathogens spreading on a network. Phys. Rev. Lett. 95 (10), 108701. Newman, M.E.J., 2006. Modularity and community structure in networks. Proc. Natl. Acad. Sci. USA 103, 8577–8582. Newman, M.E.J., 2009. Random graphs with clustering. Phys. Rev. Lett. 103 (5), 058701. Newman, M.E.J., Girvan, M., 2004. Finding and evaluating community structure in networks. Phys. Rev. E 69, 026113. ˜ o, M., Feng, Z., Marcheva, M., Castillo-Chavez, C., 2005. SIAM. J. Appl. Math. 65, Nun 964–982. Palla, G., Dere´nyi, I., Farkas, I., Vicsek, T., 2005. Uncovering the overlapping community structure of complex networks in nature and society. Nature 435, 814–818. Reichardt, J., Bornholdt, S., 2006. Statistical mechanics of community detection. Phys. Rev. E 74, 016110. Salathe´, M., Jones, J.H., 2010. Dynamics and control of diseases in networks with community structure. PLoS Comput. Biol. 6 (4), e1000736. Souza, S.R., Goncalves, S., A dynamical model for competing opinions, /arXiv:1201.1572v1S [physics.soc-ph]. Thadakamalla, H.P., Albert, R., Kumara, S.R.T., 2005. Search in weighted complex networks. Phys. Rev. E 72 (6), 066128. Trapman, P., 2007. On analytical approaches to epidemics on networks. Theor. Popul. Biol. 71 (2), 160–173. Trpevski, D., Tang, W.K.S., Jocarev, L., 2010. Model for rumor spreading over networks. Phys. Rev. E 81, 056102. Volz, E.M., 2011. Dynamics of infectious disease in clustered networks with arbitrary degree distributions. PLoS Comput. Biol. 7, e1002042. Volz, E.M., Miller, J.C., Galvani, A., Ancel Meyers, L., 2011. Effects of heterogeneous and clustered contact patterns on infectious disease dynamics. PLoS Comput. Biol. 7 (6), e1002042. Wang, Y., Xiao, G., Liu, J., 2012. Dynamics of competing ideas in complex social systems. New J. Phys. 14, 013015. Watts, D.J., Strogatz, S.H., 1998. Collective dynamics of ‘small-world’ networks. Nature 393 (6684), 440–442. Wilf, H.S., 1994. Generating functionology, 2nd ed. Academic Press, Inc. Zhou, J., Liu, Z., Li, B., 2007. Influence of network structure on rumor propagation. Phys. Lett. A 6, 458–463.