SIS epidemic propagation on hypergraphs

# SIS epidemic propagation on hypergraphs

Ágnes Bodó, Gyula Y. Katona, Péter L. Simon
###### Abstract

Mathematical modeling of epidemic propagation on networks is extended to hypergraphs in order to account for both the community structure and the nonlinear dependence of the infection pressure on the number of infected neighbours. The exact master equations of the propagation process are derived for an arbitrary hypergraph given by its incidence matrix. Based on these, moment closure approximation and mean-field models are introduced and compared to individual-based stochastic simulations. The simulation algorithm, developed for networks, is extended to hypergraphs. The effects of hypergraph structure and the model parameters are investigated via individual-based simulation results.

Institute of Mathematics, Eötvös Loránd University Budapest, Hungary

Numerical Analysis and Large Networks Research Group,

Department of Computer Science and Information Theory,

Budapest University of Technology and Economics, Budapest, Hungary

Keywords: SIS epidemic; mean-field model; exact master equation, hypergraph

AMS classification: 05C65, 60J28, 90B15, 92D30

corresponding author

email: simonp@cs.elte.hu

## 1 Introduction

Spreading processes on networks has a well-developed theory, based on the simple idea that the higher connectivity of the network enhance the spreading process. In the language of epidemic propagation this means that infection pressure on a susceptible individual increases with the number of infectious neighbours. Typically it is assumed that the probability of infection is proportional to the number of infectious neighbours. This idea leaded to different propagation models from exact master equations [16, 17] to different types of mean-field models, such as homogeneous and heterogeneous pairwise models [8, 9, 10], effective degree models [13], edge based compartmental models [14] and individual based models [15, 18], mentioning only the most widely used ones. It is also well-known that the community structure has strong impact on the spread of the epidemic. There are many publications studying household structure, workplaces, schools.

A real modeling framework of epidemic propagation has to take into account that (a) the community is built up from small units, such as households and workplaces and (b) the infection pressure on a susceptible individual in a unit is not proportional to the number of infected individuals. For example, the increase of the number of infected individuals from 5 to 10 in a workplace with 20 individuals does not necessarily mean that the infection probability is doubled in that place. The aim of this paper is to develop the theory of epidemic propagation on hypergraphs that enables us to model both the nonlinear dependence of the infection pressure and the community structure. The modeling paradigm is that the population consists of nodes of a hypergraph that is given by the hyperedges, which are simply subsets of the vertex set. For example, each household and each workplace can be considered as a hyperedge. The usual graphs can be considered as the special case of hypergraphs, when each hyperedge consists of two vertices that form an edge of the graph. In the widely-used propagation models the rate of infection for a susceptible node in a household is , where is the number of its infectious neighbours in the household and is the per contact infection rate. In our new model, to be developed in this paper, it is with some possibly non-linear function defined later. This function describes that the infection pressure is discounted as the number of infectious individuals is increased. If an individual belongs to two or more hyperedges, e.g. to a household and to a workplace, then the rate of infection is the sum of the rates in each hyperedge, e.g. , if the individual has infected neighbours at home and at the workplace. It is important to note that each hyperedge can be replaced by a clique (fully connected complete graph) and then the propagation can be considered on a graph in the conventional way by discounting the infection pressure for large number of infectious neighbours. In that case the infection rate will be instead of , hence it cannot be distinguished that somebody has many infectious neighbours at home and only a few in the workplace, or a moderate number at both places.

Besides developing the theory of epidemic propagation on hypergraphs we will run simulations on hypergraphs, derive the exact master equations and present a mean-field ODE approximation. For the purpose of simulations we created hypergraphs in several different ways. On one hand, the set of vertices was once partitioned into households and once into workplaces randomly, in this case each node belongs to exactly two hyperedges. As a second alternative, a Barabási-Albert random graph was generated and the cliques were determined by a suitable algorithm. These cliques were considered as the hyperedges of the newly created hypergraph. This is motivated by the fact that in data mining it is common practice to use different algorithms for identifying cliques as communities from given data of binary relations (i.e. a graph). Finally, we extended the configuration model to generate random hypergraphs with prescribed number and size of hyperedges.

The effect of community structure on the propagation process has been studied in the literature recently, since non-trivial community structure occurs not only in epidemiology but also in other systems in biology, computer science and engineering. Epidemic spreading on networks with overlapping community structure is considered in [7], allowing that the rate of infection is different in the communities. Clique networks describe the phenomenon of individuals attending different groups in the form of multiple clique types in [19], where theoretical analysis with a predefined clique degree distribution is performed. The authors of [11] conclude that modeling with hypergraphs may result in a more precise description of biological processes, moreover, they foresee that ”applications of hypergraph theory in computational biology will increase in the near future”. A voter-like model is investigated for interacting particle systems on hypergraphs in [12]. In their model large blocks of vertices may flip simultaneously, since vertices in a hyperedge change simultaneously their opinion to the majority opinion of the hyperedge.

The paper is structured as follows. In Section 2 our model for epidemic propagation on hypergraphs is introduced. Different methods for creating hypergraphs are presented in Section 3. The simulation results for epidemic propagation on different hypergraphs are shown in Section 4. In Section 5 we extend our theory [16] for deriving exact master equations to the case of propagation on hypergraphs. Exact and approximating mean-field equations for the expected value of the number of infected nodes are derived and compared to simulations in Section 6.

## 2 Model formulation

Our mathematical model starts from a hypergraph that consists of a set of nodes and a set of hyperedges where each hyperedge is a subset of , that is for all . The pair is called a hypergraph. The nodes in our model represent the individuals and the hyperedges are the units of the community structure such as households or workplaces. The methods for creating hypergraphs will be dealt with in the next section, here denotes a hypergraph in general. The process considered in this paper is (susceptible-infected-susceptible) epidemic propagation. This means that each node may be in one of the two states susceptible or infected/infectious. These states will be denoted by and , respectively. A susceptible individual can become infectious after contacting infectious ones and an infectious one can recover and become susceptible again after some time (not depending on the states of its neighbours). Both infection and recovery are governed by a Poisson processes. This means that an infected individual recovers with probability in a small time interval , where is called the recovery rate. A susceptible individual becomes infected with probability in a small time interval . The rate is given as

 r=τ∑hf(kh),

where the summation is for those hyperedges that contain the susceptible node, denotes the number of infected nodes in the hyperedge and is a given function. We note that in the well-known case of propagation on graphs is the identity function, i.e. , yielding the infection rate in the form , where is the total number of infected neighbours. The choice of function is artificial in this paper. The typical functional form is an inverse tangent like function that is close to the identity around zero and becomes constant for large values of its argument. The behaviour of the propagation process can be tested numerically for a large class of functions. A careful analysis shows that the qualitative behaviour of the process can be understood by using the simplest piece-wise linear function

 f(x)={x, if 0≤x≤cc, if x>c. (1)

This function is parametrized by a single parameter , which can be interpreted as a threshold value. If the number of infected nodes in a hyperedge is smaller than , then the infection is proportional to the number of infected neighbours as in the conventional network case, while above this threshold value the number of infected nodes in a hyperedge does not increase the infection pressure. We note that the desirable approach for quantitative analysis would be to determine the functional form and parameters of by fitting it to real epidemic propagation data. This approach is beyond the scope of this paper and may be the subject of future work.

Thus our model is specified by choosing a hypergraph and the function . Then the full mathematical model is a continuous time Markov chain with a state space of size , since to determine the state of the system it has to be given for each node if it is susceptible or infected. The hyperedges of the hypergraph and the function determine the transition probabilities of the Markov chain as it will be presented in Section 5. The full set of master equations of this Markov chain are practically impossible to solve because of the extremely large size of the system, hence individual-based simulations are carried out and the dependence of the system behaviour on the hypergraph structure and on the function is investigated via simulations, the results of which are shown in Section 4. Before turning to simulation results let us review the methods that were used to create our hypergraphs.

## 3 Hypergraph types used in the simulations

As we have mentioned in the introduction simulations were preformed on hypergraphs generated in three different ways. In this section we introduce these models.

### 3.1 Bi-uniform hypergraph model

In the first model, each individual (a vertex of ) belongs to precisely two edges of the hypergraph, one will correspond to the household of the individual and the other to the workplace of the person. It is assumed that the households are disjoint, and each person works in one workplace. Thus the set of edges corresponding to the households, , form a partition of . For sake of simplicity we also assume that each household has precisely members, thus is a -uniform hypergraph consisting of disjoint edges of size . With similar assumptions, the edge set of the workplaces, , is a partition of into -element sets. We generate the partitions and randomly: always pick the next element to the next edge with uniform probability from the remaining unpartitioned elements. Finally we take the union of the two uniform hypergraphs, and obtain a special bi-uniform hypergraph .

### 3.2 BA-cliques model

The second model is based on the preferential attachment model invented by Barabási and Albert [1]. They described a randomized algorithm that constructs graphs that are nowadays considered to be one of the best models for numerous natural and human-made structures, including social networks. There are several slightly different ways of implementing this algorithm, however in each case one obtains a graph with similar characteristics.

In [4] it was shown that the number of vertices with degree at least is , and in [5] it was proved that the diameter of such a graph is asymptotically . These results match the measures of many natural networks.

In our case, first we generated a random graph using the above method with our own implementation, then it was converted to a hypergraph on the same vertex set. The graph itself represents only the structure of binary relations, however, we need a representation of the structure of some communities, consisting of more members. We assume that each member of a community is in relation with each other member of the community, so members of a community forms a complete subgraph (i.e. there is en edge between any pair of members). Therefore to find these communities we listed all complete subgraphs of our graph using the software called CFinder [6], these are the edges of our hypergraph. In this way we obtain many small size hyperedges and few large ones [2].

It is worth mentioning that this algorithm has an exponential worst case running time for general graphs as there may be exponentially many cliques. Fortunately one can show that for graphs obtained by the preferential attachment model this is very unlikely, and the algorithm is most likely efficient in this case.

### 3.3 Configuration model

In this model our aim is to generate a random, uniform, regular hypergraph on given number of vertices. Therefore, we fix the number of vertices (), the size of each hyperedge () and the number of edges each vertex is contained in (i.e. the degree of each vertex, ). Simple double counting shows that in this way the number of hyperedges is . In a more general setup, one can start with a prescribed size for each hyperedge individually, and a prescribed degree for each vertex. For example, we can set half of the edges to have size and the other half size , while half of the vertices have degree and the other half (see Fig. 6.)

All hypergraphs can be associated with a bipartite graph (i.e. a 2-colorable graph) in the following way. Let be the set of vertices of the hypergraph, this will be one of the color classes in the bipartite graph. The other color class will consist of vertices corresponding to each edge of the hypergraph (). A vertex in is adjacent to vertex in if the vertex in contained in edge corresponding to the vertex in . Since in our first case every vertex is contained in hyperedges, each vertex in has degree , and since each edge contains vertices, every vertex in has degree . In the general case, vertices of have the prescribed degrees, and the degree of the vertices in will be the prescribed size of the corresponding edge.

Thus it is enough to generate such a random bipartite graph. For this we used the configuration model of Bollobás [3]. For each vertex in the bipartite graph, take the prescribed number of half-edges, then randomly combine two such half edges from the opposite side. In this way parallel edges may occur. It is shown in [3] that the expected number of parallel edges is relatively small, especially if the prescribed degrees are small compared to . So we simply chose to delete any parallel edges. As a result a few edges of the corresponding hypergraph contain less than vertices, and a few vertices are contained in less than edges. However, it has a negligible effect on the simulation results. The situation is similar in the general case.

## 4 Simulation results

### 4.1 Description of the simulation algorithm

Let denote the state of the system at time . Its -th coordinate, , is , if the -th node is susceptible and , if it is infected. The hypergraph is given by its incidence matrix , the rows of which correspond to the nodes and the coloumns correspond to the hyperedges. That is if node belongs to the -th hyperedge and it is zero otherwise. Then the product gives the number of infected nodes in the different hyperedges, that is its -th coordinate, , is the number of infected nodes in the -th hyperedge. Since infection and recovery are governed by Poisson processes, a susceptible individual, which has infected neighbours in hyperedge , becomes infected with probability in a small time interval . Similarly, an infected individual recovers with probability in a small time interval . Thus, assuming that node is susceptible at time , i.e. , the rate at which it becomes infected is

 τM∑j=1Jijf((x(t)J)j).

We apply the widely used individual-base stochastic simulation as for conventional networks. At a given time instant a vector is generated with random numbers, then the algorithm runs through all the nodes from to . If the -th node is susceptible, i.e. , then it becomes infected at time , if

 ri<1−exp(−τM∑j=1Jijf((x(t)J)j)Δt).

If the -th node is infected, i.e. , then it becomes susceptible at time , if

 ri<1−exp(−γΔt).

This process is run with sufficiently small time steps until the final time is reached. We note that running the simulation by using the Gillespie algorithm we obtain completely similar results, when is chosen sufficiently small. Then several simulations, started with the same initial condition, are averaged. The simulation results are dealt with in the next subsections.

### 4.2 The effect of function f

Here it is studied how the propagation process is affected by the choice of the function . As it was declared in Section 2, in this paper we use the function given in (1). This function is parametrised by a single parameter , which can be interpreted as a threshold value in the following sense. If the number of infected nodes in a hyperedge is smaller than , then the infection is proportional to the number of infected neighbours as in the conventional network case, while above this threshold value the number of infected nodes in a hyperedge do not increase the infection pressure.

Consider first the case of bi-uniform hypergraphs constructed from households and workplaces introduced in Subsection 3.1. A random hypergraph was constructed with households of size and workplaces of size , where each node belongs to exactly two hyperedges, a household and a workplace. Then simulations were run with recovery rate and infection rate for different values of . The time dependence of the number of infected nodes is shown in Figure 2 for three different values of . One can see that for greater values of we face stronger epidemic as it is expected. In order to compare these simulation results to the usual simulation on graphs we created a graph from the hypergraph by exchanging hyperedges to cliques, i.e. to complete subgraphs. In other words, each node in a hyperedge is connected to every other node of that hyperedge with usual edges. It may happen that two nodes are in the same household and in the same workplace, then two edges are created between them in the course of the above process. In this case the two edges are weighted with , in order to make the new network comparable to the original hypergraph. This way a weighted graph is created from the hypergraph. Then simulations were run on this weighted graph as well and the time dependence of the number of infected nodes is compared to that obtained from the simulations on the hypergraph. Figure 2 clearly shows that for large values of the process on the hypergraph is basically the same as the process on the corresponding conventional (weighted) network. This is explained by the simple fact that for the discount effect of the function cannot come to play. Namely, there are no more than nodes in the hyperedges, hence for those values of the number of infected neighbours that can occur in this hypergraph. As it was already mentioned in Section 2, in the course of a quantitative approach this Figure would enable us to fit the value of parameter by comparing real data to the curves obtained from simulations for different values of .

Consider now the case of hypergraphs created from the cliques of Barabási-Albert networks, see Subsection 3.2. We constructed a network with nodes by using the preferential attachment model with . Then the cliques were determined and substituted by hyperedges. Simulations were run with recovery rate and infection rate for different values of . The time dependence of the number of infected nodes is shown in Figure 3 for three values of . One can again see that for greater values of we face stronger epidemic as it is expected. The smallest value, , is below the average hyperedge size, hence the effect of the function can be clearly seen. Namely, the infection is smaller than on the weighted graph, which is created from the hypergraph by changing hyperedges to cliques. On the other hand, for the largest value, , which is greater than the average size of the hyperedges, the hypergraph structure has negligible effect in the steady state compared to the conventional propagation on the weighted network. More importantly, we can see that the hyperedge model has strong effect in the early stage of the epidemic, when large cliques are infected probably, for which the hyperedge size is larger than these values of .

Consider now the case of random hypergraphs created by the configuration model, see Subsection 3.3. We constructed a regular hypergraph, in which all hyperedges contain nodes and the degree of each node is , i.e. each node belongs to hyperedges. Then simulations were run with recovery rate and infection rate for different values of . The time dependence of the number of infected nodes is shown in Figure 4 for two values of . One can again see that for greater values of we face stronger epidemic as it is expected. The smaller value, , is below the hyperedge size, hence the effect of the function can be observed. Namely, the infection is smaller than on the weighted graph, which is created from the hypergraph by changing hyperedges to cliques. On the other hand, for the greater value , which is the same as the size of the hyperedges, the hypergraph structure has no effect compared to the conventional propagation on network, since for those values of the number of infected neighbours that can occur in this hypergraph.

### 4.3 The effect of the structure of the hypergraph

Now we fix the function , that is fix a value of and investigate how the parameters of the hypergraph affect the spreading process. Consider first the case of bi-uniform hypergraphs constructed from households and workplaces in Subsection 3.1. A random hypergraph was constructed with households of size and workplaces of size , where each node belongs to exactly two hyperedges, a household and a workplace. Then simulations were run with recovery rate and infection rate for different values of and . Figure 5 shows that the number of infected nodes during the spreading process depends in a non-trivial way on the values of and . One can observe that increasing the size of the hyperedges, for example, comparing the case , to the case , , the infection becomes stronger. On the other hand, comparing the case , to the case , we can say that on a more heterogeneous (in the sense of hyperedge sizes) hypergraph the propagation starts faster, but it can settle at a smaller steady state value than in the case of a more homogeneous hypergraph.

Let us turn to the study of the effect of degree heterogeneity. A hypergraph is constructed with hyperedges, half of them is of size and and the other half is of size . Half of the nodes has degree , i.e. belong to hyperedges, and the half of them is of degree . We note that the conservation relation must hold. In Figure 6 the time dependence of the number of infected nodes is shown for different degree distributions. We can see that for a more homogeneous degree distribution the spreading starts slower but it ends at a higher value. It is also important to note that for the most heterogeneous case, when the size of the small hyperedges is less than , the discount effect of the function applies only for the large hyperedges, while in the most homogeneous case with , the function affects all hyperedges.

The effect of heterogeneity is also shown in Figure 7, where the spread on a regular and a bimodal hypergraph is compared to that on a hypergraph with five different hyperedge sizes. The hypergraph have nodes and consist of hyperedges. The number of hyperedges is the same in each size category, i.e. for the bimodal hypergraph there are 200 hyperedges with both sizes and and for the hypergraph with five hyperedge sizes there are 80 hyperedges with sizes (). One can again observe that the fastest spread is on the most heterogeneous hypergraph in the initial phase of the process, while this network leads to the least prevalence in the final stage of the process.

## 5 Master equation

Here the exact master equations of SIS epidemic propagation on an arbitrary hypergraph with nodes are derived. The master equations form a linear system of ordinary differential equations, the coefficients of which are the transition rates from one state to another.

The state space of the process is containing elements, since each node can be in one of two states or . The state space is divided into classes according to the number of infected nodes. Let denote the state where every node is susceptible, i.e. . Let be the subset with states having infected nodes, containing states. Finally, let be the state in which every node is infected, i. e. .

The elements of are denoted by . Let be the type of the -th node in the state , so that the value of is either or . The state of the system can change in two ways:

1. Infection: a susceptible node becomes infected, this is an transition, where and are such that there exists , for which , and for every . Furthermore, node has an infected neighbour, which can be formally expressed as follows: there exists an , such that and and are in the same hyperedge.

2. Recovery: an infected node becomes susceptible, this is an transition, where and are such that there exists , for which , and for every .

Let denote the probability that the system is in state at time . Let

 Xk(t)=(Xk1(t),Xk2(t),…,Xkck(t))

denote the probability of states containing infected nodes, where . The above transitions define a system of linear differential equations with constant coefficients for , known as the Kolmogorov equations or master-equations. The number of the infectious nodes changes by one at most in each time step, thus the master-equation can be written in the following block tridiagonal form:

 ˙Xk=AkXk−1+BkXk+CkXk+1,k=0,1,…,N, (2)

where and are zero matrices. In matrix form:

 ˙X=PX,

where

 P=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝B0C00000A1B1C10000A2B2C20000A3B3C30⋮⋮⋯⋯⋯⋮00⋯⋯ANBN⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

The matrices describe infection and describe recovery. The structure of the network is reflected by the matrices .

Let denote the -th element of , which shows the transition rate from state to state . The class has elements and the class consists of terms, therefore the matrix has rows and columns. The entry is non-zero if and only if and differ only at one position. Let the -th node be this one, i.e. , and for every . Furthermore, there exists a number , such that and is in the same hyperedge as . Then the transition rate is given as

 Aki,j=τ∑h:l∈hf(Nh(Sk−1j)), (3)

where denotes the number of infected nodes in hyperedge in the state and is the function given in (1). Summing these equations for one obtains for every that

 ck∑i=1Aki,j=τNfSI(Sk−1j), (4)

where

 NfSI(Skj)=∑l:Skj(l)=S∑h:l∈hf(Nh(Skj)). (5)

That is denotes the sum of the values for susceptible nodes in state .

Let denote the element in the -th row and -th column of the matrix . This gives the transition rate from state to state . The class has elements and the class has terms, therefore has rows and columns. The entry is non-zero if and only if and differ only at one node. Let this node be the -th one, that is , and for every . In this case . The number of infected nodes in state is , hence there are elements in the -th column of the matrix , which are equal to , all other entries are zero. Thus for every we have

 ck∑i=1Cki,j=γ(k+1). (6)

The matrices are diagonal with rows and coloumns. The elements of denote the rate of the type transition, which are non-zero if and only if . Since the column-wise sum of the elements of are zero, the matrix is determined as follows:

 Bki,i=−ck+1∑j=1Ak+1j,i−ck−1∑j=1Ck−1j,i. (7)

As an example, the master equations are determined below for the hypergraph with nodes shown in Figure 1. The incidence matrix of the hypergraph is

 J=⎛⎜ ⎜ ⎜⎝100110011101⎞⎟ ⎟ ⎟⎠.

In this case the state space has elements and we have four classes according to the number of infected nodes:

 X0 =XSSSS, X1 =(XSSSI,XSSIS,XSISS,XISSS), X2 =(XSSII,XSISI,XSIIS,XISSI,XISIS,XIISS), X3 =(XSIII,XISII,XIISI,XIIIS), X4 =XIIII.

The vector of probabilities is and the matrix of the linear system of ODEs yielding the master equations is

 P=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝B0C0000A1B1C1000A2B2C2000A3B3C3000A4B4⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠.

The submatrices are be given as follows

 A1=⎛⎜ ⎜ ⎜⎝0000⎞⎟ ⎟ ⎟⎠,A2=τ⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝f(1)f(1)00f(1)0f(1)00f(1)f(1)0f(1)00f(1)000000f(1)f(1)⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,
 A3=τ⎛⎜ ⎜ ⎜ ⎜⎝2f(1)2f(1)2f(1)000f(1)00f(1)2f(1)00f(2)0f(2)0f(2)00f(1)02f(1)f(1)⎞⎟ ⎟ ⎟ ⎟⎠,
 A4=τ(f(2),f(1)+f(2),2f(1),f(1)+f(2)),
 B0=(0),C0=(γ,γ,γ,γ),
 C1=⎛⎜ ⎜ ⎜ ⎜⎝γγ0γ00γ0γ0γ00γγ00γ000γγγ⎞⎟ ⎟ ⎟ ⎟⎠,
 C2=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝γγ00γ0γ0γ00γ0γγ00γ0γ00γγ⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠,C3=⎛⎜ ⎜ ⎜ ⎜⎝γγγγ⎞⎟ ⎟ ⎟ ⎟⎠.

For example, the elements of the third row of the matrix are

 (0,τf(2),0,τf(2),0,τf(2)),

which describe the following rates of transitions: , , , , , . The first, third and fifth transitions can not be realized, as the states differ in more than one node. During the second transition, the first node becomes infected, hence we need to compute the term . This is equal to , because the first node is in a hyperedge, where it has two infected neighbours. The other transitions can be determined in a similar way.

Using this method one can determine the master equations for an arbitrary hypergraph theoretically. However, we have to note that the master equations are useful mainly from the theoretical point of view, as the above construction can only be carried out practically for hypergraphs of moderate size, because of the huge number of equations. This motivates the derivation of approximating systems, called mean-field equations that will be dealt with in the next section.

## 6 Mean-field theory

### 6.1 Exact differential equations for the expected number of infected and susceptible nodes

The main idea of mean-field theory is to consider some expected quantities instead of the probabilities of each individual state. The most important quantities are the expected number of infected and susceptible nodes that can be given as

 [I](t)=N∑k=0kck∑j=1Xkj(t),[S](t)=N∑k=0(N−k)ck∑j=1Xkj(t). (8)

During the derivation of differential equations for these expected values below, we will need the quantity

 [SI]=N∑k=0ck∑j=1NfSI(Skj)Xkj(t), (9)

that can be considered as the generalization of the average number of edges defined in conventional networks. We remind that is defined in (5). Now we are in the position to derive the exact differential equations for the expected values and starting from the master equations (2).

###### Theorem 1

The expected values and satisfy the following differential equations for an arbitrary hypergraph.

 ˙[S] =γ[I]−τ[SI], (10) ˙[I] =τ[SI]−γ[I]. (11)

Proof

Introducing the notation we have

 ck∑j=1Xkj=SkXk

hence (8) takes the form

 [I](t)=N∑k=0kSkXk,[S](t)=N∑k=0(N−k)SkXk. (12)

Equation (7) can be written as

 Bki,i=−(Sk+1Ak+1)i−(Sk−1Ck−1)i,

and using that is a diagonal matrix we get

 Bki,i=(SkBk)i.

Thus for every

 (SkBk)i=−(Sk+1Ak+1)i−(Sk−1Ck−1)i

holds that can be written as

 Sk+1Ak+1+SkBk+Sk−1Ck−1=0, (13)

holding for every , where and are zero matrices.

Differentiating the function and using the derivatives of given by (2) leads to

 ˙[I]=N∑k=0kSk˙Xk=N∑k=0kSk(AkXk−1+BkXk+CkXk+1)= =N∑k=1kSkAkXk−1+N∑k=0kSkBkXk+N−1∑k=0kSkCkXk+1= =N−1∑k=0(k+1)Sk+1Ak+1Xk+N∑k=0kSkBkXk+N∑k=1(k−1)Sk−1Ck−1Xk= =N∑k=0((k+1)Sk+1Ak+1+kSkBk+(k−1)Sk−1Ck−1)Xk.

Thus from equation (13) we obtain the following differential equation:

 ˙[I]=N∑k=0(Sk+1Ak+1−Sk−1Ck−1)Xk.

Now, the desired equation, (11) can be derived by using the proposition below. The proof for is similar.

###### Proposition 1

The matrices and satisfy the following identities.

Proof

From equation (6) one obtains

 (Sk−1Ck−1)j=ck−1∑i=1Ck−1i,j=γk,

for every , therefore , proving the first part of the statement.

The second part follows from the first one by using equation (12).

To prove the last statement write equation (4) as

 (Sk+1Ak+1)j=ck+1∑i=1Ak+1i,j=τNfSI(Skj),

for an arbitrary . Therefore

 N∑k=0Sk+1Ak+1Xk=N∑k=0ck∑j=1(Sk+1Ak+1)jXkj=τN∑k=0ck∑j=1NfSI(Skj)Xkj(t)=τ[SI]

that we wanted to prove.

Theorem 1 provides exact differential equations for the expected number of susceptible and infected nodes, however, these differential equations are not self-contained, since the function is not known. The next step of the mean-field theory is to derive an approximation of in terms of and in order to make the system closed. This is called a moment closure approximation, which is the subject of the next subsection.

### 6.2 Closed mean-field equations and their comparison to simulation

Our aim now is to derive an approximation to given in (9), which is based on the definition in (5).

Consider first hypergraphs describing the network structure with households of size and workplaces of size . In these hypergraphs, a node belongs to exactly two hyperedges, a household and a workplace. Denoting the total number of infected nodes by , the average number of infected neighbours of a node can be approximated by in a household and by in a workplace. Thus the second summation in (5) consists of two terms and is approximated as

 f(H−1NI)+f(W−1NI).

Since this is independent of , the double sum in (5) reduces to

 NfSI(Skj)≈(N−k)[f(H−1NI)+f(W−1NI)],

because the number of susceptible nodes in state is . Now (9) leads to the approximation

 [SI]≈[f(H−1NI)+f(W−1NI)]N∑k=0ck∑j=1(N−k)Xkj(t)=
 [f(H−1NI)+f(W−1NI)](N−[I]).

Thus equation (11) can be approximated as

 ˙I=τ(N−I)[f(H−1NI)+f(W−1NI)]−γI, (14)

The solution of this equation is compared to simulation in Figure 8 for two different hypergraphs. We can observe that the mean-field approximation gives better agreement when the hypergraph is homogeneous, i.e. for the case , .

Consider now regular random hypergraphs, in which each node belongs to hyperedges and each hyperedge is of size . Denoting the number of infected nodes by , the average number of infected neighbours of a node in a hyperedge can be approximated by . Thus the second summation in (5) consists of a single term and can be approximated as . Since this is independent of , the double sum in (5) reduces to

 NfSI(Skj)≈(N−k)df(e−1NI),

because the number of susceptible nodes in state is and each node belongs to hyperedges. Now (9) leads to the approximation

 [SI]≈df(e−1NI)N∑k=0ck∑j=1(N−k)Xkj(t)=df(e−1NI)(N−[I]).

Thus, for regular random hypergraphs, equation (11) can be approximated as

 ˙I=τ(N−I)df(e−1NI)−γI. (15)

The solution of this equation is compared to simulation in Figure 9 for two different values of . We can see that for regular random hypergraphs the mean-field approximation performs well and for the steady state it gives excellent agreement. Similarly to the case of networks, the solution of the mean-field equation increases faster than the simulated curve, but their steady states are close to each other.

## 7 Discussion

In this paper the methods of mathematical modeling of epidemic propagation on networks is extended to hypergraphs. The aim of this extension is to account for both the community structure and the nonlinear dependence of the infection pressure on the number of infected neighbours. The novelty of the model is that a susceptible individual is assumed to become infected with probability in a small time interval with rate , where the summation is for those hyperedges that contain the susceptible node, denotes the number of infected nodes in the hyperedge and is a function given in (1). The simulation algorithm, developed for networks, is extended to hypergraphs to account for this new transition probability function. Individual-based stochastic simulations were run on three different types of hypergraphs, configuration random hypergraphs, hypergraphs with hyperedges created from the cliques of a power-law random graph and bi-uniform hypergraphs the vertices of which are divided into households and workplaces randomly. The effects of hypergraph structure and the model parameters are investigated via individual-based simulation results. The exact master equations of the epidemic spreading are derived for an arbitrary hypergraph given by its incidence matrix. Based on these, moment closure approximation and mean-field models are introduced and compared to individual-based stochastic simulations.

The paper is the first step in extending the mathematical modeling of epidemic spreading from networks to hypergraphs. There are several directions where this extension can be continued. One of these is to develop and investigate pairwise models, in which not only the expected value of susceptible and infected nodes, but also the extension of pairs, as introduced in (9), are determined. The exact master equation enables us to use lumping for reducing the size of the system when the network has a special symmetry. The idea of lumping is also extendable to the case of hypergraphs. Our results show that the simple mean-field models developed in Section 6.2 performs well for homogeneous hypergraphs. The investigation of heterogeneous hypergraphs and the derivation and investigation of the corresponding heterogeneous mean-field models is also a challenging subject. The functional form and parameters of in (1) could be determined based on real epidemic propagation data by using some fitting procedure. This approach is beyond the scope of this paper and may be the subject of future work.

## Acknowledgements

Péter L. Simon acknowledges support from Hungarian Scientific Research Fund, OTKA, (grant no. 115926). Gyula Y. Katona acknowledges support from OTKA (grant no. 108947).

## References

• [1] R. Albert, A.L. Barabási, Statistical mechanics of complex networks, Reviews of Modern Physics 74 (1) (2002) 47–97.
• [2] G. Bianconi, M. Marsili, Emergence of large cliques in random scale-free networks, Europhysics Letter 74 (2006) 740–746.
• [3] B. Bollobás, A probabilistic proof of an asymptotic formula for the number of labelled regular graphs, European Journal of Combinatorics 1 (4) (1980) 311–316.
• [4] B. Bollobás, O.M. Riordan, J. Spencer, G. Tusnády, The degree sequence of a scale-free random graph process, Random Structures and Algorithms 18 (2001) 279–290.
• [5] B. Bollobás, O. Riordan, The Diameter of a Scale-Free Random Graph, Combinatorica 24 (1) (2004) 5–34.
• [6] G. Palla, I. Derényi, I. Farkas, T. Vicsek, Uncovering the overlapping community structure of complex networks in nature and society, Nature 435 (2005) 814–818. Software: CFinder 2.0.6, http://www.cfinder.org
• [7] J. Chen, H. Zhang, Z-H. Guan, T. Li, Epidemic spreading on networks with overlapping community structure, Physica A: Statistical Mechanics and its Applications 391 (2012) 1848–1854.
• [8] L. Danon, A.P. Ford, T. House, C.P. Jewell, M.J. Keeling, G.O. Roberts, J.V. Ross, M.C. Vernon, Networks and the Epidemiology of Infectious Disease, Interdisciplinary Perspectives on Infectious Diseases 2011 (2011) 1–28.
• [9] J.P. Gleeson, High-accuracy approximation of binary-state dynamics on networks, Physical Review Letters 107 (2011) 1–9.
• [10] T. House, M.J. Keeling, Insights from unifying modern approximations to infections on networks, Journal of The Royal Society Interface 8 (2011) 67–73.
• [11] S. Klamt, U-U. Haus, F. Theis, Hypergraphs and cellular networks, PLoS Computational Biology 5 (5) (2009) e1000385.
• [12] N. Lanchier, J. Neufer, Stochastic dynamics on hypergraphs and the spatial majority rule model, Journal of Statistical Physics 151 (2012) 21–45.
• [13] J. Lindquist, J. Ma, P. van den Driessche, F.H. Willeboordse, Effective degree network disease models, Journal of Mathematical Biology 62 (2011) 143–164.
• [14] J.C. Miller, A.C. Slim, E.M. Volz, Edge-based compartmental modelling for infectious disease spread, Journal of The Royal Society Interface 9 (70) (2012) 890–906.
• [15] K.J. Sharkey, Deterministic epidemiological models at the individual level, Journal of Mathematical Biology 57 (2008) 311–331.
• [16] P.L. Simon, M. Taylor, I.Z. Kiss, Exact epidemic models on graphs using graph automorphism driven lumping, Journal of Mathematical Biology 62 (2010) 479–508.
• [17] P. Van Mieghem, J. Omic, R. Kooij, Virus spread in networks, IEEE/ACM Transactions on Networking 17 (2009) 1–14.
• [18] P. Van Mieghem, The N-intertwined SIS epidemic network model, Computing 93 (2011) 147–169.
• [19] B. Wang, L. Cao, H. Suzuki, K. Aihara, Impacts of clustering on interacting epidemics, Journal of Theoretical Biology 304 (2012) 121–130.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters