# Contagion Processes on Time-Varying Networks with Homophily-Driven Group Interactions

Authors: Li Ding; Ping Hu
Journal: Complexity
 (2019)

---

## Abstract

The complicated interaction patterns among heterogeneous individuals have a profound impact on the contagion process in the networks. In recent years, there has been increasing evidence for the emergence of many-body interactions between two or more nodes in a wide range of biological and social networks. To encode these multinode interactions explicitly, the simplicial complex is now a popular alternative to simple networks. Meanwhile, the time-varying network has been acknowledged as a key ingredient of the contagion process. In this paper, we consider the connectivity pattern of networks affected by the homophily effect associated with individual attributes and investigate the impact of homophily-driven group interactions on the contagion process in temporal networks. The simplicial complex modeling framework is adopted to capture stochastic interactions between passively selected nodes in the paradigm of activity-driven networks. We study the evolution of infection and the epidemic threshold of the contagion process by both analytical and numerical methods. Our results on statistical topological properties of instantaneous network may shed light on accurately characterizing the evolution curve of infection. Furthermore, we show the impact of the homophily-driven interaction pattern on the epidemic threshold, which generalizes the existing results on both the paradigmatic activity-driven network and the simplicial activity-driven network.


---

## Body

## 1. Introduction

Network modeling plays a critical role in identifying structural properties and analyzing contagion processes on networks, such as the spreading of epidemic and malware, as well as the diffusion of news and ideas. Meanwhile, it has been acknowledged that the intrinsic characteristic differences among individuals and the complicated connectivity patterns are two main important factors determining the properties of the contagion process [1–3]. For instance, in a social network, individuals are distinguished by their attributes such as gender and age. These, together with biased interaction patterns, induce heterogeneous rates of adoption of an idea [4]. Consequently, there are substantial works dedicated to the comprehension of the typical behavioral characteristics and distinctive connectivity features, such as the homophily [4, 5], reinforcement effect [6, 7], heterogeneous degree distribution [8–10], and community structures [11, 12].Generally, there is a tacit hypothesis in network modeling approaches that the structure of complex systems is reducible to the pairwise interactions (links) of their entities (nodes) [13]. In recent years, the research community has accumulated overwhelming evidence for the emergence of many-body interaction patterns in a vast body of biological and social systems. For instance, the rich datasets available have revealed that interactions occur between more than two nodes in many systems [14]. These complicated interactions cannot be encoded in the simple network and thus prompt the need for characterizing generalized network structures. Simplicial complexes are a generalization of networks that describe interactions between more than two nodes. They can involve any number of nodes. For instance, simplices of dimension d=0,1,2, and 3 are nodes, links, triangles, and tetrahedra, respectively, while d-dimensional simplices are their d-dimensional generalizations. A simplicial complex is a finite collection of simplices of different dimensions properly glued together. It has recently been adopted in modeling many complex interacting systems, such as the brain [15], biological protein-interaction systems [16], and social systems [17].In many cases, the interactions among system entities are rapidly changing and evolving over time. The understanding of time-varying networks and their impacts on the dynamical process has been a long-standing area of research. Recently, an activity-driven (AD) modeling framework has been put forward to address the time-scale separation between the evolving network structure and the dynamical process [18]. In this framework, each node is characterized by a specific activity rate. At each time step, a node may become active and generate finite connections with other nodes randomly selected from the network. The simple formulation of AD networks and its extensions is amenable to analytical treatment and thus facilitates the study of contagion processes on time-varying networks. Examples include epidemic spreading processes with different network structures [19, 20], information diffusion processes with complex adoption models [21, 22], and the interplay between the epidemic spreading and awareness diffusion with different risk perception mechanisms [23–25].However, the above studies focused on the pairwise interactions between individuals, while the evolving simplicial interactions are rarely considered. This gap is bridged by the simplicial activity-driven (SAD) model proposed in the recent work [26]. Within this framework, a multiagent interaction is depicted by a simplex of nodes, which implies the interaction between any two nodes in a group created by an active node. Although this simplex description is accurate for some systems such as scientific collaboration networks, a more general group interaction may contain different interaction patterns among individuals. Take the microblogging network as an example, when a user actively posts a text and interacts with its followers, two or more of its followers may also be connected because they share common interests. Thus, the potential role of possible interactions among nodes connected to a common active node in the epidemic spreading process deserves further discussions.As mentioned above, the intrinsic characteristic of an individual is critical for shaping the individual interaction pattern. Inspired by the existence of temporal homophily revealed in the temporal motifs of communication networks [27], we consider the case in which the interaction patterns are affected by the individual attribute. We assume that there are two types of nodes in the network depending on their attributes. At each step, an active node randomly chooses its neighbors, and these neighbor nodes would interact at different probabilities determined by their types. In particular, nodes tend to interact with others of the same type. That is to say, the formation of group interactions is homophily-driven. This model extends the paradigms of typical AD models and the SAD model, which either ignore any interaction among neighbor nodes or assume complete interactions among neighbor nodes. We refer to this model as the complex activity-driven (CAD) model since the group interactions facilitated by active nodes are now described by a simplicial complex.In this paper, the effects of homophily-driven group interactions on the contagion processes in time-varying networks are investigated. We consider the susceptible-infected-susceptible (SIS) model for the epidemic spreading process and provide analytical analysis for both the evolution process and the epidemic threshold. In majority of existing studies, the characterization of both the evolution and the steady state of infection relies on extensive Monte Carlo simulations. Once the epidemics outbreaks, the spreading models based on the classical heterogeneous mean-field (HMF) method can only provide an upper bound for the infection level. Here, we propose a new model by developing the classical HMF method to more accurately predict the evolution of the epidemic curve. To this end, we explore the topological properties of the instantaneous network generated at each step and obtain the degree distribution and the connection correlations among nodes. It is shown that introducing these statistical information of the instantaneous network contributes to accurately characterize the evolution of epidemic. Then, we derive the epidemic threshold with the analytical statistical model. The threshold is related to the connectivity pattern determined by the node’s activity potential and the homophily effect, as well as the proportion of two types of nodes. Our analytical results can cover and generalize the existing results in both classical AD and SAD networks.The rest of this paper is constructed as follows. In Section2, the basic CAD model is established and the topological properties of the instantaneous network are analyzed. In Section 3, the statistical SIS epidemic model is formulated and the epidemic threshold is discussed. In Section 4, numerical simulations for the epidemic spreading process are presented. Conclusions are drawn in Section 5.


## 2. Model

This section introduces the CAD model and investigates the topological properties of the instantaneous CAD network.


### 2.1. Complex Activity-Driven Network

The complex activity-driven network is an extension of the paradigmatic activity-driven network that incorporates group interactions at each step. In the framework of AD models, each node (individual)i in a population of N nodes is endowed with an activity rate ai, which is taken from a predefined activity probability distribution Fa. ai is defined as the probability per unit time of node i to create interactions with other nodes. In addition to the heterogeneous activity rates, individuals in the population are distinguished by their attributes. We assume that there are two types of nodes, A and B. The proportions of type A nodes and type B nodes in the network are p and 1−p, respectively. For simplicity, we assume that the activity rates of nodes are irrelevant to their types, and thus the activity probability distributions for both type A and B nodes are expected to be Fa. Similar to the approach for static networks, we group nodes into different classes according to their types and activities. Nodes with the same type and the same activity belong to the same class and are considered to be statistically equivalent. We denote the class NaA as the collection of all nodes of type A with the activity rate a. Other classes are denoted similarly.The instantaneous networkGt of the CAD model is generated as follows: At each time step t, Gt starts with N disconnected nodes. With probability aiΔt, a node i becomes active and generates connection with m other nodes chosen randomly (as shown in Figure 1). All the selected nodes are neighbors of i, which together with the node i, form the group i. There are two different roles in this group: center and neighbor, which are played by the active node i and all other nodes connected by i, respectively. The neighbor nodes may interact with each other, and the probabilities for an interaction between two neighbor nodes with the same type and different types are λ1 and λ2, respectively. Since the interactions between neighbor nodes are stochastic, there may be different interaction patterns in one group, such as pairwise interactions and interactions involving more than two nodes. With one center node and m neighbor nodes in each group, the dimension of group interactions may vary from 1 to m. This is different from the SAD model that considers an m-simplex in each group. The instantaneous network Gt is then modeled as a simplicial complex composed of lower dimension simplicial complexes in each group. All these simplicial complexes are deleted at the next time step t+Δt. Without loss of generality, we set Δt=1 in the following.Figure 1
The CAD model. In the network, there are two types of nodes, with the blue (green) circles representing typeA (B). At each time step, a type A node i becomes active (labeled with solid circle) with the probability aiΔt. Once the node i is active, it interacts with m neighbor nodes randomly chosen from all the other nodes in the network. The purple (orange) dotted line denotes the possible interaction between two neighbor nodes with the same type (different types). With the given group interaction probabilities, a simplicial complex with one link and two triangles is generated at this step.


### 2.2. Topological Properties of the Instantaneous Network

To quantitatively characterize the epidemic spreading process, we analyze the degree distribution and the connection correlation for nodes belonging to different classes in the instantaneous network. SinceGt evolves with time, the topological properties of Gt are analyzed in a statistical way. Based on the definition of m, ai, and Fa, it is expected that there are Na active nodes with Na groups in Gt, where a=∑aaFa. In each group, there are at most mm−1/2 edges among m neighbor nodes. Since these m neighbor nodes are chosen randomly by the center node, an interaction between any two neighbor nodes occurs with the average probability λ1p2+1−p2+2λ2p1−p. Thus, it is expected that there are Eg edges in each group, with(1)Eg=m1+λ1p2+1−p2+2λ2p1−p2m−1.From the perspective of directed networks, the number of edges existing inGt is expected to be E=2NaEg.Next, we discuss the degree distribution of nodes belonging to classNaA in the instantaneous network. At each step, the inactive nodes may be passively connected by the active nodes. The probability PA,ainactk that an inactive node i∈NaA has degree k is(2)PA,ainactk=∑k1=k/mkPak1Pgk−k1,where x denotes the minimum integer that is no less than x, Pak1 denotes the probability that the node i is connected by k1 active nodes in Gt, and Pgk−k1 denotes the probability of k−k1 interactions between the node i and other neighbor nodes. The expression of Pak1 is the same as that for AD networks [24], which is given by(3)Pak1=maNk11Nk11−1NN−k1.In the thermodynamic limit, i.e.,N⟶∞, Pak1 can be rewritten as(4)Pak1=mak1k1!e−ma.In addition, the expression ofPgk−k1 is(5)Pgk−k1=k1m−1k−k1λAk−k11−λAk1m−k,with λA=λ1p+λ21−p denoting the average probability that the node i of type A interacts with a neighbor node in the group.In equation (2),k/m implies that the inactive node i with degree k must be connected by at least k/m active nodes. This is due to the fact that a node can have at most m connections in one group. Once the node i is connected by k1≥k/m active nodes, there are k1 edges between the node i and the k1 center nodes in the network. Accordingly, the rest of k−k1 edges of the node i are created by group interactions with neighbor nodes. Since the proportions of the nodes of type A and B are, respectively, p and 1−p, the average probability for an interaction between the node i and a neighbor node in the group is expected to be λA. Thus, the k−k1 interactions between node i and all other k1m−1 neighbor nodes included in k1 groups occur with the probability Pgk−k1.For each active node with degreek, there must be m edges emitted from it. Hence, other k−m edges are generated by other active nodes or through group interactions. Similar to the inactive case, the probability PA,aactk that an active node i∈NaA has degree k is(6)PA,aactk=0,k<m,∑k1=k−m/mkPak1Pgk−m−k1,k≥m.Combining equations (2) and (6), the instantaneous degree distribution PaAk of nodes in class NaA is given by(7)PaAk=aPA,aactk+1−aPA,ainactk.The analysis for the nodes of typeB is analogous to that of type A. For the nodes of type B with activity rate a in class NaB, we can easily derive the instantaneous degree distributions PaBk as(8)PaBk=aPB,aactk+1−aPB,ainactk,with λA in PA,aactk and PA,ainactk replaced by λB=λ2p+λ11−p in PB,aactk and PB,ainactk.In Figure2, we compare the statistical results (7) and (8) with those obtained by Monte Carlo simulations. Due to the heterogeneous activity rates in the real world, we consider different activity rates among individuals in the network. The activity probability distribution is chosen as Fa∝a−2. For simplicity, we assume that there are n different activity rates in the network and they are uniformly sampled in 0,1 [24]. The sampling interval is 1/N, and the minimum activity rate in the network is denoted as a1. For the example of n=5, the activity rate is selected from a∈0.2,0.4,0.6,0.8,1 with a1=0.2. In this way, the heterogeneity of node activity rates is adjustable by the value of n. The degree distribution of the nodes with activity rate a1 is taken as an illustration example and presented in Figure 2. It is shown that the theoretical results of both Pa1Ak and Pa1Bk are consistent with the simulation results. To further measure the error between the theoretical and simulation results, we define an evaluation parameter, i.e.,(9)σn=∑k∑i=1nPaiAk−P^aiAk+PaiBk−P^aiBk,with PaiAk and PaiBk from equations (7) and (8) and their counterparts P^aiAk and P^aiBk obtained by simulations. According to the value of σn in Figure 2, the absolute error between theoretical and simulation results is very small. This error is reasonable due to the finite network size and independent realizations.Figure 2
The degree distribution of the nodes with activitya1 in the instantaneous network Gt. Node types in (a) and (b) are, respectively, A and B. The dash-dotted line (circle), solid line (square), and dotted line (diamond) represent the theoretical (simulation) results for n=5, n=10, and n=20, respectively. For the three candidates of n, the values of σn are 0.0572, 0.1736, and 0.5779, respectively. Other parameters are set to be N=10,000,p=0.3,m=5,λ1=0.9,andλ2=0.1. The simulation results are averaged over 1,000 realizations.
(a)
(b)In addition to the degree distribution, we need to analyze the connection correlation among different node classes in the instantaneous networkGt. The connection correlation describes a preference for nodes to interact with others in terms of their class. It is similar to the degree-degree correlation in the static network [8] and the activity-activity correlation in AD networks [24]. Here, we refer to it as the class-class correlation since both the type and the activity rate of a node are considered in each class.To this end, we defineFNaA,Na′A as the probability of finding a node in class NaA and another node in class Na′A at two ends of an edge randomly chosen in Gt. Also, FnNaA is defined as the probability that there is a node belonging to the class NaA at the end of a randomly selected edge in Gt. Obviously, FNaA,Na′A represents the ratio of the number of edges ENaA,Na′A created by one node in class NaA and the other node in Na′A, to the total number of edges E in the instantaneous network Gt. According to the generating rules of Gt, we can obtain(10)ENaA,Na′A=pNFamapFa′+pNFa′ma′pFa+NampFam−1pFa′λ1=p2NmFaFa′a+a′+λ1am−1.The first (second) term corresponds to an active node in classaa′ that becomes active and connects to a node in class a′a. The third term describes the case where both a node in class a and a node in class a′ are chosen as neighbor nodes of an active node. Similarly, the number of edges ENaA,Na′B related to classes NaA and Na′B can be written as(11)ENaA,Na′B=p1−pNmFaFa′a+a′+λ2am−1,with equations (1), (10), and (11), and the explicit expressions of FNaA,Na′A=ENaA,Na′A/E and FNaA,Na′B=ENaA,Na′B/E can be easily obtained. The analysis of FNaB,Na′A and FNaB,Na′B is carried out in a similar way. Obviously, FNa′A,NaA=FNaA,Na′A, FNa′B,NaA=FNaB,Na′A=FNaA,Na′B=FNa′A,NaB, FNa′B,NaB=FNaB,Na′B,∀a,a′. In addition, the sum rule of ∑a,a′FNaA,Na′A+2FNaA,Na′B+FNaB,Na′B=1 is satisfied.According to the definition ofFnNaA, FnNaA equals the ratio of the number of edges ENaA connected to the nodes in class NaA to the total number of edges E in the instantaneous network Gt. That is,(12)ENaA=pNFama+pNFama1+m−1λ1p+λ21−p=pNmFaa+a+m−1aλ1p+λ21−p.The first term corresponds to an active node in classNaA that generates m edges connected to m other randomly selected nodes. The second term stems from the fact that an inactive node i in class NaA is connected by an active node in Gt, and as the role of neighbor in the group, node i interacts with other neighbor nodes with certain probabilities. Moreover, the explicit expression of FnNaA can be obtained with equations (1) and (12). Based on equations (10)–(12), it is easily to prove that the quantities FNaA,Na′A, FNaA,Na′B, and FnNaA are satisfied with FnNaA=∑a′FNaA,Na′A+FNaA,Na′B.We defineFNa′ANaA as the conditional probability that a randomly selected edge of the nodes in class NaA points to a node in the class Na′A. Then, combining equations (10)–(12), FNa′ANaA can be written as(13)FNa′ANaA=FNaA,Na′AFnNaA=pFa′a+a′+λ1am−1a+a+m−1aλ1p+λ21−p,which indicates a connection correlation between the class NaA and Na′A in the instantaneous network Gt. Following the analysis steps of equation (13), we can derive the following conditional probabilities for other three cases:(14)FNa′BNaA=1−pFa′a+a′+λ2am−1a+a+m−1aλ1p+λ21−p,FNa′ANaB=pFa′a+a′+λ2am−1a+a+m−1aλ2p+λ11−p,FNa′BNaB=1−pFa′a+a′+λ1am−1a+a+m−1aλ2p+λ11−p.The combination of any two activity ratesa and a′ can be analyzed by taking similar steps (10)–(14):In Figure3, the conditional probabilities of nodes with the minimum activity rate a1 are presented. Here, FNaA/BNa1A FNaA/BNa1B denotes the conditional probability with which an edge of nodes in the class Na1A (Na1B) is connected to a node in the class NaA or NaB. It can be found that the simulation results match well with the theoretical results (13) and (14). Also, it is observed that the homophily-driven group interactions play a significant role in shaping the connection pattern of nodes. Particularly, as shown in Figure 3(a), nodes in class Na1A tend to interact with nodes in class Na1A rather than nodes in class Na1B, even though the latter is dominant in quantity. Similarly, we introduce an evaluation parameter εn to further measure the error between theoretical and simulation results. That is,(15)εn=∑i=1n∑j=1nΔFNajANaiA+ΔFNajBNaiA+ΔFNajANaiB+ΔFNajBNaiB,with ΔFNajCNaiD=FNajCNaiD−F^NajCNaiD,C,D=A,B and FNajCNaiD from equations (13) and (14) and its counterpart F^NajCNaiD obtained by simulations. According to Figure 3, εn is small under the three candidates of n.Figure 3
The conditional probability corresponding to the nodes with activitya1 in the instantaneous network Gt. In (a) and (b), we show FNaA/BNa1A and FNaA/BNa1B for the type A nodes and type B nodes, respectively. Theoretical and simulation results are labeled in the same manner as in Figure 2. For n=5,10, and 20, εn=0.015, 0.0702, and 0.3981, respectively. Other parameters are set to be the same as those in Figure 2.
(a)
(b)


## 2.1. Complex Activity-Driven Network

The complex activity-driven network is an extension of the paradigmatic activity-driven network that incorporates group interactions at each step. In the framework of AD models, each node (individual)i in a population of N nodes is endowed with an activity rate ai, which is taken from a predefined activity probability distribution Fa. ai is defined as the probability per unit time of node i to create interactions with other nodes. In addition to the heterogeneous activity rates, individuals in the population are distinguished by their attributes. We assume that there are two types of nodes, A and B. The proportions of type A nodes and type B nodes in the network are p and 1−p, respectively. For simplicity, we assume that the activity rates of nodes are irrelevant to their types, and thus the activity probability distributions for both type A and B nodes are expected to be Fa. Similar to the approach for static networks, we group nodes into different classes according to their types and activities. Nodes with the same type and the same activity belong to the same class and are considered to be statistically equivalent. We denote the class NaA as the collection of all nodes of type A with the activity rate a. Other classes are denoted similarly.The instantaneous networkGt of the CAD model is generated as follows: At each time step t, Gt starts with N disconnected nodes. With probability aiΔt, a node i becomes active and generates connection with m other nodes chosen randomly (as shown in Figure 1). All the selected nodes are neighbors of i, which together with the node i, form the group i. There are two different roles in this group: center and neighbor, which are played by the active node i and all other nodes connected by i, respectively. The neighbor nodes may interact with each other, and the probabilities for an interaction between two neighbor nodes with the same type and different types are λ1 and λ2, respectively. Since the interactions between neighbor nodes are stochastic, there may be different interaction patterns in one group, such as pairwise interactions and interactions involving more than two nodes. With one center node and m neighbor nodes in each group, the dimension of group interactions may vary from 1 to m. This is different from the SAD model that considers an m-simplex in each group. The instantaneous network Gt is then modeled as a simplicial complex composed of lower dimension simplicial complexes in each group. All these simplicial complexes are deleted at the next time step t+Δt. Without loss of generality, we set Δt=1 in the following.Figure 1
The CAD model. In the network, there are two types of nodes, with the blue (green) circles representing typeA (B). At each time step, a type A node i becomes active (labeled with solid circle) with the probability aiΔt. Once the node i is active, it interacts with m neighbor nodes randomly chosen from all the other nodes in the network. The purple (orange) dotted line denotes the possible interaction between two neighbor nodes with the same type (different types). With the given group interaction probabilities, a simplicial complex with one link and two triangles is generated at this step.


## 2.2. Topological Properties of the Instantaneous Network

To quantitatively characterize the epidemic spreading process, we analyze the degree distribution and the connection correlation for nodes belonging to different classes in the instantaneous network. SinceGt evolves with time, the topological properties of Gt are analyzed in a statistical way. Based on the definition of m, ai, and Fa, it is expected that there are Na active nodes with Na groups in Gt, where a=∑aaFa. In each group, there are at most mm−1/2 edges among m neighbor nodes. Since these m neighbor nodes are chosen randomly by the center node, an interaction between any two neighbor nodes occurs with the average probability λ1p2+1−p2+2λ2p1−p. Thus, it is expected that there are Eg edges in each group, with(1)Eg=m1+λ1p2+1−p2+2λ2p1−p2m−1.From the perspective of directed networks, the number of edges existing inGt is expected to be E=2NaEg.Next, we discuss the degree distribution of nodes belonging to classNaA in the instantaneous network. At each step, the inactive nodes may be passively connected by the active nodes. The probability PA,ainactk that an inactive node i∈NaA has degree k is(2)PA,ainactk=∑k1=k/mkPak1Pgk−k1,where x denotes the minimum integer that is no less than x, Pak1 denotes the probability that the node i is connected by k1 active nodes in Gt, and Pgk−k1 denotes the probability of k−k1 interactions between the node i and other neighbor nodes. The expression of Pak1 is the same as that for AD networks [24], which is given by(3)Pak1=maNk11Nk11−1NN−k1.In the thermodynamic limit, i.e.,N⟶∞, Pak1 can be rewritten as(4)Pak1=mak1k1!e−ma.In addition, the expression ofPgk−k1 is(5)Pgk−k1=k1m−1k−k1λAk−k11−λAk1m−k,with λA=λ1p+λ21−p denoting the average probability that the node i of type A interacts with a neighbor node in the group.In equation (2),k/m implies that the inactive node i with degree k must be connected by at least k/m active nodes. This is due to the fact that a node can have at most m connections in one group. Once the node i is connected by k1≥k/m active nodes, there are k1 edges between the node i and the k1 center nodes in the network. Accordingly, the rest of k−k1 edges of the node i are created by group interactions with neighbor nodes. Since the proportions of the nodes of type A and B are, respectively, p and 1−p, the average probability for an interaction between the node i and a neighbor node in the group is expected to be λA. Thus, the k−k1 interactions between node i and all other k1m−1 neighbor nodes included in k1 groups occur with the probability Pgk−k1.For each active node with degreek, there must be m edges emitted from it. Hence, other k−m edges are generated by other active nodes or through group interactions. Similar to the inactive case, the probability PA,aactk that an active node i∈NaA has degree k is(6)PA,aactk=0,k<m,∑k1=k−m/mkPak1Pgk−m−k1,k≥m.Combining equations (2) and (6), the instantaneous degree distribution PaAk of nodes in class NaA is given by(7)PaAk=aPA,aactk+1−aPA,ainactk.The analysis for the nodes of typeB is analogous to that of type A. For the nodes of type B with activity rate a in class NaB, we can easily derive the instantaneous degree distributions PaBk as(8)PaBk=aPB,aactk+1−aPB,ainactk,with λA in PA,aactk and PA,ainactk replaced by λB=λ2p+λ11−p in PB,aactk and PB,ainactk.In Figure2, we compare the statistical results (7) and (8) with those obtained by Monte Carlo simulations. Due to the heterogeneous activity rates in the real world, we consider different activity rates among individuals in the network. The activity probability distribution is chosen as Fa∝a−2. For simplicity, we assume that there are n different activity rates in the network and they are uniformly sampled in 0,1 [24]. The sampling interval is 1/N, and the minimum activity rate in the network is denoted as a1. For the example of n=5, the activity rate is selected from a∈0.2,0.4,0.6,0.8,1 with a1=0.2. In this way, the heterogeneity of node activity rates is adjustable by the value of n. The degree distribution of the nodes with activity rate a1 is taken as an illustration example and presented in Figure 2. It is shown that the theoretical results of both Pa1Ak and Pa1Bk are consistent with the simulation results. To further measure the error between the theoretical and simulation results, we define an evaluation parameter, i.e.,(9)σn=∑k∑i=1nPaiAk−P^aiAk+PaiBk−P^aiBk,with PaiAk and PaiBk from equations (7) and (8) and their counterparts P^aiAk and P^aiBk obtained by simulations. According to the value of σn in Figure 2, the absolute error between theoretical and simulation results is very small. This error is reasonable due to the finite network size and independent realizations.Figure 2
The degree distribution of the nodes with activitya1 in the instantaneous network Gt. Node types in (a) and (b) are, respectively, A and B. The dash-dotted line (circle), solid line (square), and dotted line (diamond) represent the theoretical (simulation) results for n=5, n=10, and n=20, respectively. For the three candidates of n, the values of σn are 0.0572, 0.1736, and 0.5779, respectively. Other parameters are set to be N=10,000,p=0.3,m=5,λ1=0.9,andλ2=0.1. The simulation results are averaged over 1,000 realizations.
(a)
(b)In addition to the degree distribution, we need to analyze the connection correlation among different node classes in the instantaneous networkGt. The connection correlation describes a preference for nodes to interact with others in terms of their class. It is similar to the degree-degree correlation in the static network [8] and the activity-activity correlation in AD networks [24]. Here, we refer to it as the class-class correlation since both the type and the activity rate of a node are considered in each class.To this end, we defineFNaA,Na′A as the probability of finding a node in class NaA and another node in class Na′A at two ends of an edge randomly chosen in Gt. Also, FnNaA is defined as the probability that there is a node belonging to the class NaA at the end of a randomly selected edge in Gt. Obviously, FNaA,Na′A represents the ratio of the number of edges ENaA,Na′A created by one node in class NaA and the other node in Na′A, to the total number of edges E in the instantaneous network Gt. According to the generating rules of Gt, we can obtain(10)ENaA,Na′A=pNFamapFa′+pNFa′ma′pFa+NampFam−1pFa′λ1=p2NmFaFa′a+a′+λ1am−1.The first (second) term corresponds to an active node in classaa′ that becomes active and connects to a node in class a′a. The third term describes the case where both a node in class a and a node in class a′ are chosen as neighbor nodes of an active node. Similarly, the number of edges ENaA,Na′B related to classes NaA and Na′B can be written as(11)ENaA,Na′B=p1−pNmFaFa′a+a′+λ2am−1,with equations (1), (10), and (11), and the explicit expressions of FNaA,Na′A=ENaA,Na′A/E and FNaA,Na′B=ENaA,Na′B/E can be easily obtained. The analysis of FNaB,Na′A and FNaB,Na′B is carried out in a similar way. Obviously, FNa′A,NaA=FNaA,Na′A, FNa′B,NaA=FNaB,Na′A=FNaA,Na′B=FNa′A,NaB, FNa′B,NaB=FNaB,Na′B,∀a,a′. In addition, the sum rule of ∑a,a′FNaA,Na′A+2FNaA,Na′B+FNaB,Na′B=1 is satisfied.According to the definition ofFnNaA, FnNaA equals the ratio of the number of edges ENaA connected to the nodes in class NaA to the total number of edges E in the instantaneous network Gt. That is,(12)ENaA=pNFama+pNFama1+m−1λ1p+λ21−p=pNmFaa+a+m−1aλ1p+λ21−p.The first term corresponds to an active node in classNaA that generates m edges connected to m other randomly selected nodes. The second term stems from the fact that an inactive node i in class NaA is connected by an active node in Gt, and as the role of neighbor in the group, node i interacts with other neighbor nodes with certain probabilities. Moreover, the explicit expression of FnNaA can be obtained with equations (1) and (12). Based on equations (10)–(12), it is easily to prove that the quantities FNaA,Na′A, FNaA,Na′B, and FnNaA are satisfied with FnNaA=∑a′FNaA,Na′A+FNaA,Na′B.We defineFNa′ANaA as the conditional probability that a randomly selected edge of the nodes in class NaA points to a node in the class Na′A. Then, combining equations (10)–(12), FNa′ANaA can be written as(13)FNa′ANaA=FNaA,Na′AFnNaA=pFa′a+a′+λ1am−1a+a+m−1aλ1p+λ21−p,which indicates a connection correlation between the class NaA and Na′A in the instantaneous network Gt. Following the analysis steps of equation (13), we can derive the following conditional probabilities for other three cases:(14)FNa′BNaA=1−pFa′a+a′+λ2am−1a+a+m−1aλ1p+λ21−p,FNa′ANaB=pFa′a+a′+λ2am−1a+a+m−1aλ2p+λ11−p,FNa′BNaB=1−pFa′a+a′+λ1am−1a+a+m−1aλ2p+λ11−p.The combination of any two activity ratesa and a′ can be analyzed by taking similar steps (10)–(14):In Figure3, the conditional probabilities of nodes with the minimum activity rate a1 are presented. Here, FNaA/BNa1A FNaA/BNa1B denotes the conditional probability with which an edge of nodes in the class Na1A (Na1B) is connected to a node in the class NaA or NaB. It can be found that the simulation results match well with the theoretical results (13) and (14). Also, it is observed that the homophily-driven group interactions play a significant role in shaping the connection pattern of nodes. Particularly, as shown in Figure 3(a), nodes in class Na1A tend to interact with nodes in class Na1A rather than nodes in class Na1B, even though the latter is dominant in quantity. Similarly, we introduce an evaluation parameter εn to further measure the error between theoretical and simulation results. That is,(15)εn=∑i=1n∑j=1nΔFNajANaiA+ΔFNajBNaiA+ΔFNajANaiB+ΔFNajBNaiB,with ΔFNajCNaiD=FNajCNaiD−F^NajCNaiD,C,D=A,B and FNajCNaiD from equations (13) and (14) and its counterpart F^NajCNaiD obtained by simulations. According to Figure 3, εn is small under the three candidates of n.Figure 3
The conditional probability corresponding to the nodes with activitya1 in the instantaneous network Gt. In (a) and (b), we show FNaA/BNa1A and FNaA/BNa1B for the type A nodes and type B nodes, respectively. Theoretical and simulation results are labeled in the same manner as in Figure 2. For n=5,10, and 20, εn=0.015, 0.0702, and 0.3981, respectively. Other parameters are set to be the same as those in Figure 2.
(a)
(b)


## 3. Analysis

In this section, we characterize the paradigmatic SIS model with the statistical properties of the instantaneous network obtained in Section2. Based on the statistical SIS model, we discuss the epidemic threshold in the CAD networks.


### 3.1. Statistical SIS Model

In the SIS model, each node can be either susceptibleS or infectious I. A susceptible node in contact with an infectious node becomes infected with an infection rate β, while each infectious node recovers from infection with a recovery rate μ.At the time stept, we denote the fractions of susceptible and infectious nodes in the class NaA as SA,at and IA,at, respectively. Then, it gives SA,at+IA,at=1. To analyze the contagion process quantitatively, we define ωA,at as the probability that an edge of nodes in class NaA points to an infectious node in Gt. Obviously, ωA,at is dependent on both the class-class connection correlation property of Gt and the infection level in each node class. Based on the analysis in Section 2, we have(16)ωA,at=∑a′FNa′ANaAIA,a′t+FNa′BNaAIB,a′t,with FNa′ANaA and FNa′BNaA given by equations (13) and (14), respectively.For a nodei in the class NaA, the probability of node i with degree k is given by the degree distribution PaAk. Based on equation (7), the fractions of susceptible and infectious nodes in the class NaA at next time step t+1 are given by(17)SA,at+1=μIA,at−SA,at∑k1−1−βωA,atkPaAk,IA,at+1=1−μIA,at+SA,at∑k1−1−βωA,atkPaAk.It is noted that the statistical properties of the instantaneous networkGt, PaAk and FNa′ANaA, are included in (17). Since PaAk and FNa′ANaA are related to the parameters p, λ1 and λ2, the potential impacts of homophily-driven group interactions on the contagion process are reflected in the statistical properties of Gt. We refer to the SIS model described by equation (17) as our statistical model. This model is different from models derived by the classical HMF theory in most of existing studies within the framework of AD models. Leveraging the HMF theory, the SIS model on the CAD networks is derived as(18)I^A,at+1=1−μI^A,at+βS^A,atma∫I^A,a′t+I^B,a′tNda′+βS^A,at∫I^A,a′t+I^B,a′tNma′da′+βS^A,at∫ma′da′×∫λ1I^A,a″t+λ2I^B,a″tNm−1da″,where S^A,at and I^A,at, respectively, denote the number of susceptible and infectious nodes in class NaA at the time step t. With S^A,at=pNFaSA,at and I^A,at=pNFaIA,at, equation (18) can be rewritten as(19)IA,at+1=1−μIA,at+βSA,atmapIAt+1−pIBt+βSA,atmpθAt+1−pθBt+βSA,atmam−1λ1pIAt+λ21−pIBt,with IAt=∑aFaIA,at, IBt=∑aFaIB,at, θAt=∑aaFaIA,at, and θBt=∑aaFaIB,at. To distinguish with the statistical model, the model (19) is referred to as the classical MF model. This model is introduced as a comparison model to measure the performance of the proposed statistical model that contains more topological properties of the instantaneous network.


### 3.2. Epidemic Threshold

To obtain the epidemic threshold, we analyze the stability aroundIA,at⟶0. In this case, we have the approximation 1−1−βωA,atk≈kβωA,at. Then, it gives(20)∑k1−1−βωA,atkPAak=∑kkβωA,atPAak=ENaApNFaβωA,at.By introducing equations (13) and (14) into ωA,at, the fraction of infectious nodes IA,at+1 in equation (17) can be rewritten as(21)IA,at+1=1−μIA,at+βSA,atmpIAta+λ1am−1+1−pIBta+λ2am−1+pθAt+1−pθBt,which is exactly equation (19) derived by the classical HMF theory. This fact then implies that the classical MF model (19) is an approximation of the statistical model (17). Meanwhile, the result of approximation is accurate only when the infection level is very low. This will be further discussed in the following numerical simulations.By integrating equation (21) over all the node classes and ignoring the second-terms in IAt, we can obtain the expression of IAt+1:(22)IAt+1=1−μIAt+βmpIAta1+λ1m−1+1−pIBta1+λ2m−1+pθAt+1−pθBt.Moreover, the closed expression forθAt+1 is derived by multiplying both sides of equation (21) by a and integrating over all classes spectrum. That is,(23)θAt+1=1−μθAt+pβmIAta2+λ1a2m−1+1−pβmIBta2+λ2a2m−1+βmapθAt+1−pθBt.Similarly, the expressions ofIBt+1 and θBt+1 that correspond to the nodes of type B can be written as(24)IBt+1=1−μIBt+βmpIAta1+λ2m−1+1−pIBta1+λ1m−1+pθAt+1−pθBt,θBt+1=1−μθBt+pβmIAta2+λ2a2m−1+1−pβmIBta2+λ1a2m−1+βmapθAt+1−pθBt.Thus, we obtain the following closed system of the master equations forIAt, IBt, θAt, and θBt:(25)∂tIAt=−μIAt+pβma1+λ1m−1IAt+1−pβma1+λ2m−1IBt+βmpθAt+1−pθBt,∂tIBt=−μIBt+pβma1+λ2m−1IAt+1−pβma1+λ1m−1IBt+βmpθAt+1−pθBt,∂tθAt=−μθAt+pβma2+λ1a2m−1IAt+1−pβma2+λ2a2m−1IBt+βmapθAt+1−pθBt,∂tθBt=−μθBt+pβma2+λ2a2m−1IAt+1−pβma2+λ1a2m−1IBt+βmapθAt+1−pθBt,whose Jacobian matrix is given by(26)J=−μ+pβma1+λ1m−11−pβma1+λ2m−1pβm1−pβmpβma1+λ2m−1−μ+1−pβma1+λ1m−1pβm1−pβmpβma2+λ1a2m−11−pβma2+λ2a2m−1−μ+pβma1−pβmapβma2+λ2a2m−11−pβma2+λ1a2m−1pβma−μ+1−pβma.The critical condition for the epidemic threshold is determined by the largest eigenvalueΛmax of J. That is, if(27)ΛmaxJ>0,the epidemics will outbreak in the network. Here, we cannot get the explicit expression of equation (27) for general systems. In the following, we will discuss two special cases of the system.Ifλ1=λ2=λ, there is no homophily effect on interactions between neighbor nodes in any groups. Thus, all nodes in the network can be considered to be of the same type. In this case, the master equations of the system can be rewritten as(28)∂tIAt=−μIAt+βma1+λm−1pIAt+1−pIBt+βmpθAt+1−pθBt,∂tIBt=−μIBt+βma1+λm−1pIAt+1−pIBt+βmpθAt+1−pθBt,∂tθAt=−μθAt+βma2+λa2m−1pIAt+1−pIBt+βmapθAt+1−pθBt,∂tθBt=−μθBt+βma2+λa2m−1pIAt+1−pIBt+βmapθAt+1−pθBt,which can be further reduced to be(29)∂tIt=−μIt+βma1+λm−1It+βmθt,∂tθt=−μθt+βma2+λa2m−1It+βmaθt,with It=pIAt+1−pIBt and θt=pθAt+1−pθBt. The corresponding Jacobian matrix J is given by(30)J=−μ+βma1+λm−1βmβma2+λa2m−1−μ+βam,and the two eigenvalues of J are(31)Λ1,2=−μ+βmλm−1+2a±βmλm−1+22a2+4a2−a22.Thus, the critical condition (27) for the presence of an endemic state is(32)βμ>2mλm−1+2a+mλm−1+22a2+4a2−a2.Particularly, withλ=λ1, the results for the cases of p=0 and p=1 are also included in (32).Furthermore, ifλ=0, no interaction between neighbor nodes will be allowed and equation (32) will recover the result β/μ>1/ma+a2 for AD networks. In addition, if λ=1, an interaction between any two neighbor nodes can be expected in each group. Then, equation (32) gives β/μ>2/mm+1a+mm+12a2+4a2−a2, which is degenerated to the results of the SAD model in [26].Next, we discuss the case where the number of nodes of both types is the same. Withp=0.5, the master equation (25) can be reduced as(33)∂tIt=−μIt+βma1+λ1+λ22m−1It+βmθt,∂tθt=−μθt+βma2+λ1+λ22a2m−1It+βmaθt.Following the analysis steps of equation (29), we can derive the eigenvalues of J determined by system (33). That is,(34)Λ1,2=−μ+βmλ1+λ2/2m−1+2a±βmλ1+λ2/2m−1+22a2+4a2−a22.Thus, the critical condition (27) is given by(35)βμ>2mλ1+λ2/2m−1+2a+mλ1+λ2/2m−1+22a2+4a2−a2.Comparing equations (32) and (35), we can see that the case of a uniform mixing of two types of nodes is equivalent to the case of λ=λ1+λ2/2 without homophily effect.


## 3.1. Statistical SIS Model

In the SIS model, each node can be either susceptibleS or infectious I. A susceptible node in contact with an infectious node becomes infected with an infection rate β, while each infectious node recovers from infection with a recovery rate μ.At the time stept, we denote the fractions of susceptible and infectious nodes in the class NaA as SA,at and IA,at, respectively. Then, it gives SA,at+IA,at=1. To analyze the contagion process quantitatively, we define ωA,at as the probability that an edge of nodes in class NaA points to an infectious node in Gt. Obviously, ωA,at is dependent on both the class-class connection correlation property of Gt and the infection level in each node class. Based on the analysis in Section 2, we have(16)ωA,at=∑a′FNa′ANaAIA,a′t+FNa′BNaAIB,a′t,with FNa′ANaA and FNa′BNaA given by equations (13) and (14), respectively.For a nodei in the class NaA, the probability of node i with degree k is given by the degree distribution PaAk. Based on equation (7), the fractions of susceptible and infectious nodes in the class NaA at next time step t+1 are given by(17)SA,at+1=μIA,at−SA,at∑k1−1−βωA,atkPaAk,IA,at+1=1−μIA,at+SA,at∑k1−1−βωA,atkPaAk.It is noted that the statistical properties of the instantaneous networkGt, PaAk and FNa′ANaA, are included in (17). Since PaAk and FNa′ANaA are related to the parameters p, λ1 and λ2, the potential impacts of homophily-driven group interactions on the contagion process are reflected in the statistical properties of Gt. We refer to the SIS model described by equation (17) as our statistical model. This model is different from models derived by the classical HMF theory in most of existing studies within the framework of AD models. Leveraging the HMF theory, the SIS model on the CAD networks is derived as(18)I^A,at+1=1−μI^A,at+βS^A,atma∫I^A,a′t+I^B,a′tNda′+βS^A,at∫I^A,a′t+I^B,a′tNma′da′+βS^A,at∫ma′da′×∫λ1I^A,a″t+λ2I^B,a″tNm−1da″,where S^A,at and I^A,at, respectively, denote the number of susceptible and infectious nodes in class NaA at the time step t. With S^A,at=pNFaSA,at and I^A,at=pNFaIA,at, equation (18) can be rewritten as(19)IA,at+1=1−μIA,at+βSA,atmapIAt+1−pIBt+βSA,atmpθAt+1−pθBt+βSA,atmam−1λ1pIAt+λ21−pIBt,with IAt=∑aFaIA,at, IBt=∑aFaIB,at, θAt=∑aaFaIA,at, and θBt=∑aaFaIB,at. To distinguish with the statistical model, the model (19) is referred to as the classical MF model. This model is introduced as a comparison model to measure the performance of the proposed statistical model that contains more topological properties of the instantaneous network.


## 3.2. Epidemic Threshold

To obtain the epidemic threshold, we analyze the stability aroundIA,at⟶0. In this case, we have the approximation 1−1−βωA,atk≈kβωA,at. Then, it gives(20)∑k1−1−βωA,atkPAak=∑kkβωA,atPAak=ENaApNFaβωA,at.By introducing equations (13) and (14) into ωA,at, the fraction of infectious nodes IA,at+1 in equation (17) can be rewritten as(21)IA,at+1=1−μIA,at+βSA,atmpIAta+λ1am−1+1−pIBta+λ2am−1+pθAt+1−pθBt,which is exactly equation (19) derived by the classical HMF theory. This fact then implies that the classical MF model (19) is an approximation of the statistical model (17). Meanwhile, the result of approximation is accurate only when the infection level is very low. This will be further discussed in the following numerical simulations.By integrating equation (21) over all the node classes and ignoring the second-terms in IAt, we can obtain the expression of IAt+1:(22)IAt+1=1−μIAt+βmpIAta1+λ1m−1+1−pIBta1+λ2m−1+pθAt+1−pθBt.Moreover, the closed expression forθAt+1 is derived by multiplying both sides of equation (21) by a and integrating over all classes spectrum. That is,(23)θAt+1=1−μθAt+pβmIAta2+λ1a2m−1+1−pβmIBta2+λ2a2m−1+βmapθAt+1−pθBt.Similarly, the expressions ofIBt+1 and θBt+1 that correspond to the nodes of type B can be written as(24)IBt+1=1−μIBt+βmpIAta1+λ2m−1+1−pIBta1+λ1m−1+pθAt+1−pθBt,θBt+1=1−μθBt+pβmIAta2+λ2a2m−1+1−pβmIBta2+λ1a2m−1+βmapθAt+1−pθBt.Thus, we obtain the following closed system of the master equations forIAt, IBt, θAt, and θBt:(25)∂tIAt=−μIAt+pβma1+λ1m−1IAt+1−pβma1+λ2m−1IBt+βmpθAt+1−pθBt,∂tIBt=−μIBt+pβma1+λ2m−1IAt+1−pβma1+λ1m−1IBt+βmpθAt+1−pθBt,∂tθAt=−μθAt+pβma2+λ1a2m−1IAt+1−pβma2+λ2a2m−1IBt+βmapθAt+1−pθBt,∂tθBt=−μθBt+pβma2+λ2a2m−1IAt+1−pβma2+λ1a2m−1IBt+βmapθAt+1−pθBt,whose Jacobian matrix is given by(26)J=−μ+pβma1+λ1m−11−pβma1+λ2m−1pβm1−pβmpβma1+λ2m−1−μ+1−pβma1+λ1m−1pβm1−pβmpβma2+λ1a2m−11−pβma2+λ2a2m−1−μ+pβma1−pβmapβma2+λ2a2m−11−pβma2+λ1a2m−1pβma−μ+1−pβma.The critical condition for the epidemic threshold is determined by the largest eigenvalueΛmax of J. That is, if(27)ΛmaxJ>0,the epidemics will outbreak in the network. Here, we cannot get the explicit expression of equation (27) for general systems. In the following, we will discuss two special cases of the system.Ifλ1=λ2=λ, there is no homophily effect on interactions between neighbor nodes in any groups. Thus, all nodes in the network can be considered to be of the same type. In this case, the master equations of the system can be rewritten as(28)∂tIAt=−μIAt+βma1+λm−1pIAt+1−pIBt+βmpθAt+1−pθBt,∂tIBt=−μIBt+βma1+λm−1pIAt+1−pIBt+βmpθAt+1−pθBt,∂tθAt=−μθAt+βma2+λa2m−1pIAt+1−pIBt+βmapθAt+1−pθBt,∂tθBt=−μθBt+βma2+λa2m−1pIAt+1−pIBt+βmapθAt+1−pθBt,which can be further reduced to be(29)∂tIt=−μIt+βma1+λm−1It+βmθt,∂tθt=−μθt+βma2+λa2m−1It+βmaθt,with It=pIAt+1−pIBt and θt=pθAt+1−pθBt. The corresponding Jacobian matrix J is given by(30)J=−μ+βma1+λm−1βmβma2+λa2m−1−μ+βam,and the two eigenvalues of J are(31)Λ1,2=−μ+βmλm−1+2a±βmλm−1+22a2+4a2−a22.Thus, the critical condition (27) for the presence of an endemic state is(32)βμ>2mλm−1+2a+mλm−1+22a2+4a2−a2.Particularly, withλ=λ1, the results for the cases of p=0 and p=1 are also included in (32).Furthermore, ifλ=0, no interaction between neighbor nodes will be allowed and equation (32) will recover the result β/μ>1/ma+a2 for AD networks. In addition, if λ=1, an interaction between any two neighbor nodes can be expected in each group. Then, equation (32) gives β/μ>2/mm+1a+mm+12a2+4a2−a2, which is degenerated to the results of the SAD model in [26].Next, we discuss the case where the number of nodes of both types is the same. Withp=0.5, the master equation (25) can be reduced as(33)∂tIt=−μIt+βma1+λ1+λ22m−1It+βmθt,∂tθt=−μθt+βma2+λ1+λ22a2m−1It+βmaθt.Following the analysis steps of equation (29), we can derive the eigenvalues of J determined by system (33). That is,(34)Λ1,2=−μ+βmλ1+λ2/2m−1+2a±βmλ1+λ2/2m−1+22a2+4a2−a22.Thus, the critical condition (27) is given by(35)βμ>2mλ1+λ2/2m−1+2a+mλ1+λ2/2m−1+22a2+4a2−a2.Comparing equations (32) and (35), we can see that the case of a uniform mixing of two types of nodes is equivalent to the case of λ=λ1+λ2/2 without homophily effect.


## 4. Numerical Simulations

In this section, we perform numerical simulations to investigate the epidemic spreading process on CAD networks. For simplicity, we consider five activity rates of the nodes, i.e.,a∈0.5,0.6,0.7,0.8,0.9 [24], and the activity rate distribution obeys Fa∝a−2. Accordingly, the proportions of nodes with the activity rate 0.5, 0.6, 0.7, 0.8, and 0.9 are 0.344, 0.239, 0.176, 0.135, and 0.106, respectively. At each time step, the number of edges generated by an active node is chosen to be m=5. In the following simulations, the size of network, the recovery rate, and the initial fraction of randomly chosen infectious nodes in the population are set to be N=5,000, μ=0.5, and I0=0.1, respectively. If not specified, we assume the proportion of type A nodes, the interaction probabilities for two neighbor nodes of the same type, and different types to be p=0.6,λ1=0.5,andλ2=0.1, respectively. All the Monte Carlo simulation results are averaged over 40 independent realizations.Firstly, we compare the performance of the statistical model and the classical MF model on characterizing the SIS epidemic spreading process. To this end, we depict the time evolution of the fraction of infectious nodes in the network. Under different infection rateβ, the evolution curves based on the statistical model (17), the classical MF model (19) and the Monte Carlo simulations are recorded. The results are presented in Figure 4. It is observed that the statistical model captures the evolution curve for all four values of β, while the classical MF model approximates the simulation results only for the case of β=0.01. This is consistent with the analysis in Section 2, as the approximation adopted in the classical MF model is accurate only when the effective infection rate is below the epidemic threshold. Furthermore, as shown in Figures 4(b)–4(d), the classical MF model provides an upper bound for the infection level in simulations. The difference between this upper bound and the actual infection level increases as the infection rate β increases. According to the results of Figure 4, we can conclude that analyzing the statistical properties of the instantaneous network is beneficial to accurately characterize the contagion process in the time-varying networks.Figure 4
The time evolution of the fraction of infectious nodes of the typeA (marked with blue), type B (marked with green), and both two types (marked with red). The circles, solid lines, and dotted lines represent the results of Monte Carlo simulations, the classical MF model, and the statistical model, respectively. (a) β = 0.01. (b) β = 0.05. (c) β = 0.1. (d) β = 0.15.
(a)
(b)
(c)
(d)Next, we investigate the impact of homophily-driven group interactions on the final steady state of the epidemic spreading. The growth of the final fraction of infectious nodesI∞ versus the infection rate β is presented in Figure 5. It is shown that I∞ decreases with the increasing p from 0 to 0.5. Since p describes the proportion of type A nodes in the network, an increasing p≤0.5 indicates more uniform mixing of both type nodes. This implies that any two neighbor nodes in any groups are more likely to be of different types. Due to the homophily effect, it is expected that there are less group interactions among nodes. Consequently, as shown in Figure 5, the contagion process in the network is impeded and the epidemic threshold is enhanced with an increasing p. Furthermore, noting that the definition of type A and type B is exchangeable, the effect of a varying p from 0.5 to 1 is the same as a varying 1−p from 0 to 0.5. This fact, together with the observation in Figure 5, demonstrates that the impact of homophily effect on the contagion process is the greatest when p=0.5. In fact, this result can be obtained by analyzing the average probability of an interaction between any two neighbor nodes. That is, for the given λ1 and λ2, the value of λ1p2+1−p2+2λ2p1−p reaches its minimum with p=0.5. In addition, as illustrated in Figure 5, the epidemic thresholds suggested by the Monte Carlo simulations and the statistical model are close to the theoretical values given by equation (27). This observation, in turn, validates our analytical results presented in Section 3.Figure 5
The final fraction of infectious nodes as a function of the infection rateβ for different proportions of type A nodes, p. The epidemic thresholds suggested by the Monte Carlo simulation and the statistical model are recorded as soon as I∞ is slightly than 0.01 (the critical value above which the realization is considered as endemic). The analytical values of the epidemic threshold given by equation (27) are labeled by yellow pentagons.Furthermore, we consider the contagion processes on the AD networks, the SAD networks, and the proposed CAD networks. According to the analysis in Section3, the AD model and the SAD model are indeed two special cases of the CAD model and they are obtained by letting λ1=λ2=0 and λ1=λ2=1, respectively.In Figure6, we depict the growth of the final fraction of infectious nodes I∞ versus the infection rate β in the three networks. Obviously, the group interactions between two neighbor nodes promote the contagion process in the network. Compared with the AD network, the epidemic threshold decreases in both CAD and SAD networks. Specifically, the critical infection rates corresponding to the epidemic threshold are satisfied with βcAD>βcCAD>βcSAD. In addition, the theoretical results of the epidemic threshold given by equations (32) and (35) are approximate to the simulation results.Figure 6
The final fraction of infectious nodes as a function of the infection rateβ in AD networks, CAD networks, and SAD networks. Here, we chose p=0.5 for the CAD networks to illustrate the theoretical results of equation (35). The epidemic thresholds suggested by the Monte Carlo simulations, the statistical model, and the theoretical analysis are recorded and labeled following the same step as in Figure 5.Finally, we investigate the impacts of the infection rateβ and the group interaction probabilities λ1andλ2 on the epidemic threshold. Figure 7 presents a comparison between the theoretical and numerical values of the epidemic threshold for different parameter combinations, where the numerically obtained steady-state infection density is color coded and shown as the background. Here, we assume λ2 is proportional to λ1, i.e, λ2=τλ1, and the proportionality factors τ are chosen as 0.1, 0.5, and 1 in Figures 7(a)–7(c), respectively. For a fixed parameter set, the dashed curve represents the theoretical value of λ1c as a function of the infection rate β based on equation (27). The markers are the corresponding numerical results of λ1c. From Figure 7, the theoretical results agree with the simulation results for all three cases of τ. For each τ, there are two critical values of β, i.e., βc1 and βc2. When β<βc1 or β>βc2, the epidemic threshold is not related to λ1. Specifically, the epidemics will vanish from the network with β<βc1 or it breaks out with β>βc2. When βc1<β<βc2, it is observed that λ1c decreases with an increase of β. Here, βc1 is obtained under λ1=1 and is determined by the value of λ2. Given the three candidates of τ in Figure 7, we can see that βc1 decreases with an increase of τ due to more group interactions between nodes with different types. Moreover, since the CAD network turns into the AD network with λ1=λ2=0, βc2 is exactly the epidemic threshold in the AD network.Figure 7
The final fraction of infectious nodes as a function of the infection rateβ and the group interaction probabilities λ1 and λ2. In each subgraph, the curve with markers represents the dependence of the predicted critical value λ1c for the onset of the epidemic on the infection rate β. λ1c suggested by Monte Carlo simulations is obtained in the same way as deriving the critical value of β in Figure 5. (a) λ2 = 0.1λ1; (b) λ2 = 0.5λ1; (c) λ2 = λ1.
(a)
(b)
(c)


## 5. Conclusions

In this paper, we investigate the effect of group interactions involving more than two individuals on the contagion process in the time-varying networks. We consider that the interactions are affected by the individual attribute. The homophily effects on interaction patterns are studied. A new network model that extends the paradigmatic activity-driven model to the framework of simplicial complex networks is presented. The statistical properties of the instantaneous network based on this model are explored and incorporated in characterizing the epidemic spreading process. Results show that these properties provide a more accurate description of the evolution of epidemic. Furthermore, it is demonstrated that the homophily-driven group interactions have a profound impact on the epidemic threshold. Our results generalize the existing results of the epidemic threshold on two paradigms of activity-driven network models. However, the proposed model also comes with some limitations. Here, we consider that the activity rate of a node is independent of its type. Future work may explore a general activity probability distribution that incorporates the effect of different node types. Moreover, we assume that the infection rates through pairwise interaction and multinode interaction are the same. This model can be further developed by exploring different infection rates for different interaction patterns.



---

*Source: 7130468-2019-10-30.xml*