Social clustering in epidemic spread on coevolving networks.

Even though transitivity is a central structural feature of social networks, its influence on epidemic spread on coevolving networks has remained relatively unexplored. Here we introduce and study an adaptive susceptible-infected-susceptible (SIS) epidemic model wherein the infection and network coevolve with nontrivial probability to close triangles during edge rewiring, leading to substantial reinforcement of network transitivity. This model provides an opportunity to study the role of transitivity in altering the SIS dynamics on a coevolving network. Using numerical simulations and approximate master equations (AMEs), we identify and examine a rich set of dynamical features in the model. In many cases, AMEs including transitivity reinforcement provide accurate predictions of stationary-state disease prevalence and network degree distributions. Furthermore, for some parameter settings, the AMEs accurately trace the temporal evolution of the system. We show that higher transitivity reinforcement in the model leads to lower levels of infective individuals in the population, when closing a triangle is the dominant rewiring mechanism. These methods and results may be useful in developing ideas and modeling strategies for controlling SIS-type epidemics.

INTRODUCTION

In recent years, the study of dynamical processes on complex networks has received significant attention in the mathematical modeling of epidemics [1,2]. There exist three distinct approaches to modeling epidemic spread as a dynamical process on networks. In the first, each node is allowed to change its state with no evolution of the underlying network structure through time, while in the second the network structure coevolves with the state of the nodes [3,4]. A third and more complex approach involves combining the inherent (disease-independent) evolution of social networks with the epidemic dynamics [5-8]. The second approach is more realistic for modeling epidemic spread in cases where the state of being infected influences the presence of contact over which the infection can spread. Pragmatically, the details of the network evolution incorporate potentially several social processes. It is a common observation that healthy individuals avoid contact with individuals suffering from an infectious disease. A recent study shows that humans can identify sick individuals through olfactory-visual cues; furthermore, the human immune system appears to use these cues to motivate healthy individuals to avoid contact with infective individuals [9], reducing the risk of contagions and increasing the biological fitness. Another mechanism that leads to the severing of ties between healthy and infective individuals is the practice of quarantine, often enforced by public health agencies to stop the spread of communicable diseases. Such a preference in attachment between individuals changes the social network and thus influences the spread of the infection. Studying this interplay between the infection spread and the coevolving network structure has the potential to provide insights into the processes of disease spread.

The features of an underlying social network can govern the propagation of an infectious disease [10-14], and there is extensive literature exploring the impact of network structure on the dynamics of epidemic spread on static networks [10-18]. For instance, Ref. [16] shows that clustering in networks could raise the epidemic threshold and degree correlations can alter the epidemic size. However, results obtained using static network models ignore effects due to the underlying social networks evolving in response to the spread of a disease. Coevolutionary network systems, also called adaptive networks, model the dynamics of node states (e.g., susceptible vs infected) together with rules for rewiring network edges in response to the observed states, providing a more general setting to model disease dynamics. Numerous studies have explored different aspects of the interplay between coevolving networks and epidemic spread [19,20], including extensions to signed networks [21]; nevertheless, a variety of open problems and challenges still exist despite increased understanding of coevolving systems [22-24].

A common problem across many coevolving network systems is that they often employ edge rewiring rules that ignore transitivity, a fundamental structural property of most networks, wherein friends of friends are more likely to be friends. For example, the typical random-rewiring rule would randomize any given network and destroy most closed triangles present in the initial conditions; this rule is very unlikely to create new triangles in the network and yields a trivial transitivity level corresponding to that expected under edge independence assumptions. Some recent attempts have been made to explore the influence of transitivity in coevolving voter models [25-27]. However, the effect of transitivity on coevolving network epidemics remains unstudied, despite the known importance of transitivity in the static network setting.

In this paper we introduce an adaptive susceptible-infected-susceptible (SIS) model on a coevolving network which incorporates a rewiring rule that reinforces transitivity similar to that employed in [26,27]. The most straightforward strategy to deal with an infectious disease is to sever all the links between susceptible and infected individuals. However, humans are social animals and need a certain amount of social relationships to function. Therefore, only discarding the infected-susceptible links without compensation is not pragmatic in all settings. Under the mechanism considered here, susceptible individuals break links to their infected neighbors, rewiring those links to either a random individual in the network or a neighbor's neighbor. The latter rewiring will lead to reinforcement of transitivity (increased clustering coefficient) in the system. While transitivity is a critical property of social networks, studies of epidemic dynamics on coevolving networks have mostly neglected this property. Our model thus provides an opportunity to study the role of transitivity in epidemic dynamics in coevolving network systems. This rewiring rule includes a parameter that governs the preferential rewiring to close triangles. We study the impact of this transitivity reinforcement on the spread of the infection and on the details of the network using numerical simulations of the process as well as an approximate master equation framework (see [28,29]). In Sec. II we describe our model and in Sec. III we present results obtained through numerical simulations, highlighting the dynamical features observed in the model. In Sec. IV we present a derivation of our semianalytical method, confirming results obtained from numerical simulations. Sections III and IV also include identification and analysis of bifurcations in the system dynamics. We conclude in Sec. V with further discussion.

MODEL

We employ the classical setting of the SIS model [30,31]. We consider a network with nodes and undirected edges (both constant), with , where is the average degree. At a given time, each node is either in a susceptible () or infected () state. Infected individuals infect each of their susceptible neighbors at rate while recovering (change to susceptible) at rate . Meanwhile, each susceptible individual breaks each edge it has with an infected neighbor at rate , rewiring this link to a susceptible node.

Generalizing previously studied models, the new link created during rewiring is made with probability to a susceptible node at distance 2, that is, a neighbor of a neighbor to whom the given node is not already connected. Otherwise (i.e., with probability ) the rewiring is made to another susceptible node selected uniformly at random from the network outside its immediate neighborhood [that is, self-loops and repeated links (multiedges) are prohibited]. Under the latter situation, uniform rewiring, there remains nonzero probability of randomly picking a neighbor of a neighbor, but this probability becomes vanishingly small for larger networks. In the event that there is no available susceptible node to rewire to (e.g., there are no distance-2 susceptible nodes), the rewiring attempt fails and the original link to the infected neighbor remains. Importantly, rewiring to a neighbor's neighbor closes a triangle between the three nodes, directly reinforcing transitivity, whereas rewiring uniformly at random tends to decrease transitivity.

NUMERICAL SIMULATION

Unless stated otherwise, our numerical results of the adaptive SIS system presented here are for networks with nodes and edges (that is, mean degree ), and 1000 simulations of the process are produced for each configuration of parameters. We performed Monte Carlo simulations with the help of largenet, a c++ library that has been widely used for simulations of large adaptive networks [32,33]. We used largenet with fixed step size , asynchronously updating one randomly selected node at a time.

Following Ref. [34], we consider the following three initial degree distributions: the Poisson distribution approximating the distribution of an Erdős-Rényi model for large ; a truncated power law distribution with and in order to make the mean degree (we also study the effect of different cutoff degrees in the Appendix A 1); and the degree regular distribution where is the Kronecker delta indicating every node has the same initial degree . We observe that at large enough time , the system approaches a statistically stationary state with features of interest such as the clustering coefficient and disease prevalence fluctuating around mean values. For the purposes of identifying this state numerically, we determine the stationary state to be reached when the level of disease prevalence (expressed as a fraction of nodes) between two consecutive integer times differs by at most . Under this definition, we find that in practice networks with the above three initial degree distributions reach their stationary states before for . In contrast, when we do not observe convergence to stationary states, even after a long time; rather, we appear to have a long, progressively slowing decay of the disease prevalence (as we will see below in Fig. 8).

FIG. 8.

(a) Disease prevalence at time . The initial degree distribution is Poisson. Different values of are indicated by color. The other parameters of the system are , and . Circles correspond to means computed over 1000 simulations. Lines are AME results. (b) Observe that for the simulations exhibit slow convergence and the AME fails to capture the correct timescale for this convergence. (c) Plot of the corresponding clustering coefficient and transitivity at for different .

To highlight the main differences between our model and previous work without reinforced transitivity [34], we first confirm that our model leads to nontrivial transitivity. We explore two closely related measures of transitivity here. The first measure is the average local clustering coefficient , where the local clustering coefficient of vertex is , with and the number of closed triplets (triangles) and the total number of triplets centered at vertex , respectively [35]. The second measure is the transitivity (often called global clustering coefficient) [36], defined as , where and are, respectively, the total number of closed triangles and the number of triplets in the network. Figures 1(a) and 1(c) show the temporal evolution of the average local clustering coefficient and the transitivity in our model, with varying from 0 to 1 in steps of 0.2, keeping all other parameters fixed and choosing the initial degree distribution to be Poisson. When , there is no reinforcement of transitivity, i.e., () for all time and we are in the same regime studied in Ref. [34]. For nonzero , the clustering coefficients rapidly increase and then slowly converge to nonzero values, except in the case , where we do not observe this convergence, and (as well as ) appears to continue to grow with a diminishing rate. We treat as a special case in our model and will study this case separately later. We also explored the influence of initial degree distributions on the final clustering coefficient. In Figs. 1(b) and 1(d) we observe that networks with the regular degree distribution initially achieve higher clustering coefficients compared to networks with the other two initial distributions. Moreover, networks with initial Poisson and truncated power law degree distributions show indistinguishable final clustering coefficients for . Finally, in Fig. 2 we examine the local clustering coefficient grouped by susceptible nodes and infected nodes at time for networks with Poisson degree distributions initially and with different values of . As increases, susceptible nodes typically achieve higher values of local clustering coefficient compared to the infected nodes.

FIG. 1.

The level of transitivity varies with time and the probability of rewiring steps made to susceptible neighbors of neighbors. The (a) average local clustering coefficient and (c) transitivity are shown versus time for simulations with initially Poisson degree distribution with mean degree . The other parameters of the system are , and . Also shown are (b) and (d) at versus for networks with the same mean degree but different initial degree distributions (Poisson), triangle; (truncated power law), circle; and (degree regular), square]. Plotted data correspond to mean values computed over 30 Monte Carlo simulations. In (a) and (c), all cases except appear to reach stationarity by time .

FIG. 2.

Local clustering coefficients grouped by susceptible nodes and infected nodes at time for (a) , (b) , (c) , and (d) . The initial degree distribution is Poisson. The other parameters of the system are , and . At time , all cases except for reach their stationary states. Each box plot is obtained from a single Monte Carlo simulation of network size . The notch of the box indicates the median (mostly zero), the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively, whiskers extend to the most extreme data points not considered as outliers, and the outliers are plotted individually (the +'s).

Given that our model demonstrates nontrivial transitivity, we now explore how the reinforced transitivity influences the evolution of disease prevalence . In Fig. 3 we plot the observed disease prevalence at long times in the parameter space, fixing the infection rate and recovery rate as in Ref. [34]. The phase diagram in Fig. 3 indicates that for large enough the system eventually becomes disease-free, presumably because susceptible nodes can easily rid themselves of infected neighbors by rewiring these edges, thus effectively quarantining the infection. When , there is no rewiring process, the network does not change, and hence has no effect. A more interesting behavior occurs for , with large values of disease prevalence that depend on in a complex manner. Since this parameter regime is most interesting for our present study, in the rest of our simulations we fix and vary to investigate the influence of transitivity reinforcement.

FIG. 3.

Phase diagrams in parameter space for the observed disease prevalence at time on networks with initial degree distributions that are (a) Poisson, (b) truncated power law, and (c) degree regular. We set here , and . At the selected time all cases except appear to have converged to their stationary states. Results are averaged over 30 simulations at each set of parameters. These diagrams are generated through bilinear interpolation from results on a regular grid, leading to some apparent grid artifacts.

Another critical phenomenon that typically occurs in coevolving network models is the fission of the underlying network as the system evolves. In Fig. 4 we explore the number of connected components () and the fraction of nodes in the largest component () at for the three initial settings of the degree distribution. We fix the mean degree , and and then compute and at . Varying from 0 to 1, we observe the increase in and decrease in , indicating that higher leads to more fragmented networks in the stationary states. Interestingly, the results corresponding to the Poisson and the truncated power law degree distributions still overlap except when .

FIG. 4.

(a) Number of components and (b) largest component fraction versus at . The initial degree distributions of the networks are Poisson (P), truncated power law (TPL), and degree regular (DR) with mean degree . The other parameters in the system are , and . Every point in the figure corresponds to the mean value computed over 30 simulations. In these simulations, the mean number of components at for Poisson, truncated power law, and degree regular distributions are 666.17, 3348.27, and 6.03, respectively, and the mean largest component fraction at for Poisson, truncated power law, and degree regular are 0.7971, 0.6828, and 0.7452, respectively. All cases except reach stationarity by time .

To further explore the dynamics of the system, we plot bifurcation diagrams in Fig. 5. The rows in Fig. 5 correspond to different initial degree distributions with columns for different values of . When , for all three initial distributions there exists a small bistable region near . The sudden jump of around the critical when (the first column of plots) implies that the transition might be discontinuous when . This discontinuity appears to shrink as increases. For , the transition appears to be continuous. While we do not conclusively confirm the nature of the transitions here, we will investigate the transitions in more detail in Sec. V. We emphasize the separation of figures corresponding to , with the axis labeled as and not , highlighting the fact that for the system does not reach a stationary state in the observed time.

FIG. 5.

Bifurcation diagrams of the disease prevalence in the stationary state versus at different values of . Each row of plots corresponds to a specific degree distribution: The first row is Poisson, the second row is truncated power law, and the third row is degree regular. Each column of plots corresponds to a particular value of , as indicated at the top of the column. Every point represents a mean of 30 realizations of the system at . The last column belongs to the case, noting that the axis is labeled (cf. ) because the simulations have not reached stationarity. The other parameters of the system are and . (We also study the effect of different in initially Poisson networks; see the Appendix A 2.)

APPROXIMATE MASTER EQUATIONS

We further study this system with approximate semianalytical techniques, defining appropriate systems of ordinary differential equations and then solving those systems numerically. The moment-closure frameworks of pair approximation (PA) and approximate master equations (AMEs) have both been used effectively in similar settings [1,28,37]. In PA, the different counts of contact pairs are used to approximate the density of triplets, thus closing the system of equations [28]. In contrast, the AME framework considers the populations of nodes according to their state, degree, and the states of their immediate neighbors, approximating other quantities from these populations. The AME is an annealed mean field approximation that deterministically approximates stochastic systems. We note in particular that with the additional variables and corresponding differential equations in AMEs, the triplet counts are precisely accounted for (cf. the closure approximations necessary to obtain triplet counts in PA) and that the AME typically provides a more accurate approximation of such dynamics and can be highly stable around the critical point of the dynamics [1,28,34,37,38] (see also [28,37,39] for comparisons between the AME and PA). We note that Ref. [38] further extended this approach to classify links according to the states, numbers of neighbors, and numbers of infected neighbors at both ends of the links. This link-based approach can further improve on the accuracy of the AME method; however, it increases the system size from to equations, where is the maximum degree. As such, we restrict our attention here to the AME method.

We employ an approach similar to that used in Ref. [34], extending the equations to include the effect of transitivity reinforcement. Let and be the fraction of susceptible and infected nodes, respectively, of degree with infected neighbors at time [see the diagram in Fig. 6(a)]. Following the notation in Ref. [34], we define the zeroth-order moments of the and distributions by the first-order moments by and the second-order moments by It is worth noting that while node states and network topologies coevolve, some quantities are conserved. For example, , since the number of nodes is fixed. Similarly, , because of the conservation of edges. We then have the following ordinary differential equation (ODE) governing the time evolution of the compartment: Similarly, the ODE for the compartment is

FIG. 6.

Illustration of an node rewiring to a neighbor's neighbor. (a) Before the rewiring, the center nodes in the diagram are of class and , respectively. Consider the two nodes involved when an edge is rewired to a distance-2 neighbor: the active node doing the rewiring and the passive node receiving the newly rewired edge. (b) Diagram representing the case when the center is passively rewired by the action of a distance-2 neighbor, with the center node moving to class . (c) Diagram representing the case when the center actively rewires to a distance-2 neighbor and the center node becomes class . The other processes of the adaptive SIS system can be similarly diagramed.

To describe the AME derivation, we will focus on explaining the equation, as the derivation of corresponding terms in the equation involves similar arguments. The terms premultiplied by and in the first line of Eq. (4) account for the center node and its neighbors recovering from and becoming infected. The term premultiplied by on the first line describes the center class in state dismissing one of its state neighbors and rewiring to a random node with state . For a more detailed discussion of the derivation of these terms, we refer the reader to Ref. [34].

The second line of Eq. (4), premultiplied by , is new to the present work, accounting for a center class with state and its distance-2 neighbors with state , breaking links to one of their neighbors and rewiring to the center node (see the diagrams in Fig. 6). To better understand the derivation of this line, we note that only nodes can rewire and be rewired to in this model (whereas neighbors are only dropped by the action of nodes). Supposing the center node is of class , then it has proportion of its neighbors infected neighbors, and we describe the number of edges among them as , making assumptions about independence to close the equations with this approximation. Similarly, the same center node class has proportion of its neighbors susceptible, with the number of edges among them described by . Therefore, the rate associated with the transition to is . Following the same arguments, we find that the corresponding rate for a node leaving the class is .

The last line of Eq. (4), premultiplied by , is the effect of a center class with state rewired to a random node in the network with state . For the compartment in Eq. (5), the center class can only lose edges due to rewiring and the value of does not change the rate at which edges to nodes are rewired. Therefore, the equation lacks terms premultiplied by and . The resulting system of coupled ordinary differential equations contains equations, where is the maximum degree.

To determine the initial conditions for the above equations, initially a fraction of randomly chosen nodes is infected, independent of node degree. This gives us the set of initial conditions where is the degree distribution at . We proceed to numerically solve the AME system of equations (4) and (5). In the preceding section we presented an overview of the general behavior of the model via explicit simulations of the model. We now carry out a detailed comparison between the simulations and AME results, confirming the accuracy of the AMEs and furthering our understanding of the adaptive SIS dynamics with reinforced transitivity.

EXPLORING THE SYSTEM THROUGH SIMULATIONS AND AMES

First we plot the disease prevalence against time with different 's and different initial degree distributions in Fig. 7. Note that [Fig. 7(a)] is the case where our model reduces to that studied in Ref. [34]. Networks with different initial degree distributions converge to different stationary disease prevalences. In the Poisson and the truncated power law cases, the disease prevalence approaches a nonzero value, i.e., the system attains an endemic state. However, in the degree regular case, gradually approaches 0, i.e., the system approaches a disease-free state. The thick lines in Fig. 7 represent results derived from the AME approach, while symbols are simulation results.

FIG. 7.

Disease prevalence at time . The initial degree distributions are Poisson (diamonds) [at ; see Eq. (1)], truncated power law (squares) [at ; see Eq. (2)], and degree regular (circles) [at ; see Eq. (3)]. Markers correspond to means computed over 1000 simulations and lines are the semianalytical AME results. The parameters used are , and , with (a) , (b) , (c) , and (d) .

Most importantly, AME captures the stationary disease prevalence in all the cases shown in the figure. Second, for small , AME approximates the temporal evolution of the disease prevalence with a high level of accuracy, even in the initial transient states, and correctly predicts the direction of the changes with . However, the differences between simulations and the AME increase at larger and are particularly obvious for the truncated power law case, but the AME still captures the overall shape of the temporal evolution of disease prevalence. (To further highlight these elements, in Fig. 21 in the Appendix A 3, we provide a zoomed-in view for the transient dynamics for some of the panels of Fig. 8.)

FIG. 21.

Time evolution of the disease prevalence , the fraction of SI links , the effective branching factor , and the average number of connections that susceptible nodes share with other susceptible nodes , for (a)–(d) and (e)–(h) . Diamond markers correspond to means computed over 200 Monte Carlo simulations (error bars would be smaller than the markers here) and the solid lines are the AME predictions. The initial degree distribution is Poisson and other parameters of the system are , and .

To study the above results in greater detail, we allow to vary from 0 to 1 while keeping all other parameters fixed and focusing on the initial Poisson degree distribution. We compare simulation and AME results in Fig. 8, illustrating one of the critical features of our model: the disease prevalence decreases with increasing . That is, we observe that systems with a stronger preference for transitivity have smaller numbers of infected individuals, all else being equal. The thick lines in Fig. 8(a) for the AME results accurately capture the final disease prevalence and roughly describe its temporal evolution for . For , the AME underestimates the stationary disease prevalence but still captures qualitative features of the evolution. Remarking again that is a special case, we note that the system does not reach a well-defined stationary state in our simulations, with continuing to decrease. In Fig. 9 we further explore the temporal evolution of when while varying from 2 to 10. We identify a slow power law decay in for ; using the curve fitting in matlab®, we find that the disease prevalence for and scales as (also shown in Fig. 9), with appearing to slowly approach a disease-free state. As the mean degree increases, the observed decay with becomes even slower, to the point that we cannot claim anything about the functional form of the long-time decay from the figure; nevertheless, we provide power law fits to the late times in the figure for comparison. Although the AME prediction for the temporal evolution for is qualitatively very different, it does appear to predict a correct disease-free stationary state (see Fig. 8).

FIG. 9.

Scaling of disease prevalence with time when . As the average degree becomes greater than 4, the scaling exponent appears to decrease in magnitude, suggesting that this phenomenon may depend on the level of sparsity of the underlying network. The parameters here are , and . Markers correspond to means computed over 30 simulations.

To summarize, we plot simulation results (markers) and AME-based estimates (lines) of versus for different initial degree distributions in Fig. 10. The AME predictions and simulation results agree very well. The role of reinforced transitivity is clear: The final disease prevalence decreases with increasing transitivity reinforcement for networks starting with Poisson and power law degree distributions. (The degree regular cases are already disease-free at stationarity for these parameters.)

FIG. 10.

Disease prevalence in the stationary state versus . The 's were calculated at time . The initial degree distributions are Poisson (diamonds) [at ; see Eq. (1)], truncated power law (squares) [at ; see Eq. (2)], and degree regular (circles) [at ; see Eq. (3)]. Markers correspond to means computed over 30 simulations and the lines are AME results. The other parameters are , and .

In addition to predictions of disease prevalence, the AME approach contains and thus predicts the degree distributions of the coevolving networks. In Fig. 11 we plot the stationary joint distributions of node states and degrees for for the three different initial distributions; that is, the degree distributions of the susceptible and infected nodes are indicated separately with normalization corresponding to the respective fraction of nodes in each state. Markers and lines in Fig. 11 represent the simulations and AME results, respectively, demonstrating that the AME provides a good estimate in each case. Note that in Fig. 11(c) the fraction of infected nodes is very close to but not exactly zero. This is because we stopped the numerical experiments when there are no edges in the system; however, there are still a few edges and isolated infected nodes and as a result we still see infected nodes with different degrees in the stationary states.

FIG. 11.

Joint distributions of node states ( and ) and degrees at stationarity. The initial degree distributions are (a) Poisson, (b) truncated power law, and (c) degree regular. Markers correspond to means computed over 200 simulations. Dashed and solid lines are the AME predictions for the and states, respectively. The parameters used are , and .

To determine how impacts the degree distribution, we plot the initial and stationary degree distributions for different 's in Fig. 12. As in previous figures, markers are simulations and lines are AME results, which appear to be in good agreement. Different colors represent different values of , with the interesting result that appears to have only a small effect on the total stationary degree distribution, even though we have seen before that does affect the disease prevalence. Indeed, for both the initial Poisson and truncated power law cases, the stationary degree distribution appears to be close to Poisson with mean degree [see the thick gray line in Figs. 12(a) and 12(b)]. This observation implies that the stationary states resulting from this process have both a Poisson degree distribution and nonzero transitivity. We recall that the model includes two kinds of rewiring: at random to any susceptible node and to a distance-2 neighbor. Neither of these has direct preference to rewire to vertices with a particular degree, except that rewiring to a distance-2 neighbor cannot rewire to a singleton. That is, both rewiring mechanisms are essentially similar to random rewiring in ignoring node degree, which leads to a Poisson degree distribution. Hence, we expect the initial degree distributions to converge to Poisson-like distributions in stationary.

FIG. 12.

Degree distributions in stationary states. The initial degree distributions are (a) Poisson, (b) truncated power law, and (c) degree regular. Asterisks in each plot are the initial degree distributions. We plot results for , 0.2, 0.4, 0.6, and 0.8, as indicated by the color bar, though in many cases the results at one value of are obscured by those at other values. The other parameters are , and . We plot results at , with all systems at these reaching their stationary states by this time. (Note that we did not plot here.) Markers correspond to means computed over 200 simulations. Lines are the AME predictions for the stationary degree distributions. The thick gray line in the background is the reference Poisson distribution expected for a network with . Observe that for both the initial Poisson and truncated power law the stationary state degree distribution falls on this thick line, implying that the final degree in both these cases is Poisson (at least approximately).

Similar to Fig. 5 in Ref. [34], we plot (Fig. 13) bifurcation diagrams of stationary disease prevalence level versus infection rate with an initial Poisson degree distribution and with different mean degrees , 4, and 8. The other parameters here include and . We obtained the average stationary disease prevalence level at each value of from 30 Monte Carlo simulations for each initial infected proportion , 0.01, 0.05, 0.99, and 0.999. We call attention to the range of the infection rate plotted in Fig. 13, ranging from 0 to 0.04 in Figs. 13(a)–13(d) for , from 0 to 0.02 in Figs. 13(e)–13(h) for , and from 0 to 0.008 in Figs. 13(i)–13(l) for . As evident in the figure, our AME calculation reasonably predicts the stationary disease prevalence level , including its phase transition and the bistability near that transition.

FIG. 13.

Bifurcation diagrams of stationary disease prevalence level versus infection rate on adaptive networks with an initial Poisson degree distribution, with for the top row, for the middle row, and for the bottom row. In each row, we plot the case when , 0.2, 0.4, and 0.6. The other parameters in the system are and . Solid lines are the predictions of our semianalytic method and the markers represent the average of 30 simulations. We ran these Monte Carlo simulations for each value , 0.01, 0.05, 0.99, and 0.999. For each run, the initial transient was discarded and the prevalence at stationarity was averaged over at least time steps.

We observed that varying (the probability of closing triangles during rewiring) can lead to the transition of the system from endemic to a disease-free state with a very slow timescale when (see Fig. 9). To further investigate whether the nature of these transitions might be a finite-size effect, in Fig. 14 we plot the susceptibility , as defined in [40], for different network sizes . With the peak of the susceptibility identifying the transition, the results in the figure show that for the network size does not greatly change the transition location in , except in the smallest networks. However, we again see a qualitative difference for , making a definitive conclusion in this case more difficult.

FIG. 14.

Susceptibility versus infected rate in adaptive networks with Poisson initial degree distributions with different network size at time for various values of : (a) , (b) , (c) , (d) , (e) , and (f) . Except for the case , all networks reached their stationary disease prevalence . The other parameters of the system are , and . Lines with different colors correspond to means computed over 30 Monte Carlo simulations of different network sizes.

In Fig. 15 we consider a complementary test of possible finite-size effects: Keeping the average degree fixed (with initial Poisson degree distributions), we vary the network size from to while considering different values of from 0 to 1. Fixing the average degree, the final disease prevalence level appears to be almost constant across network sizes. Nevertheless, we note that the values considered in Fig. 15 are clearly above the transition between disease-free and endemic cases, though the used with in Fig. 15(b) is much closer to the transition than that in Fig. 15(a). Combined with Fig. 14, where the susceptibility and phase transition seem to not depend on the network size, these findings suggest that, in general, the stationary disease prevalence level does not change appreciably with the size of the system here. While we acknowledge that the precise details in and immediately around the transition in might be expected to be more sensitive to system size, such an exploration would involve an even more intensive numerical investigation simultaneously considering variations in , and . Meanwhile, the above results together appear to support the hypothesis that the qualitative nature of the transitions observed herein at various parameter settings are not likely to be consequences of finite system sizes.

FIG. 15.

Stationary disease prevalence level (as measured at stopping time ) versus network size on adaptive networks with an initial Poisson degree distribution with (a) and (b) . In (a) the infected and rewiring rates chosen are: and . In (b) the infected and rewiring rates chosen are and . The other parameters in the system are and . Circles (and error bars) represent the mean (and standard deviation) of outcomes from 30 Monte Carlo simulations. We vary from 0 to 1, recalling that does not reach stationarity by . We note that is consistent across network sizes in the region tested in this paper.

Furthermore, to demonstrate the influence of on the transitions, in Fig. 16(a) we plot a phase diagram for the observed stationary disease prevalence in the parameter space . We observe that as increases, the critical value of where the system changes from disease-free to endemic decreases. To supplement the phase diagram, in Fig. 16(b) we plot the numerically computed derivative of disease prevalence along the direction, , observing only a narrow range with a large positive value near the transition points apparent in Fig. 16(a). Also, the jump in occurs at a lower as increases, again showing the dependence of the transitions on . Finally, when , we see no jumps of when varies, further suggesting that the system behavior at might be qualitatively different. The observed transitions might also depend on other parameters in the system. In the Appendix A 2, we study the influence of (initial fraction of infected nodes) on the transitions, and the effect of remains.

FIG. 16.

(a) Phase diagram in parameter space for the observed disease prevalence at time on the networks with Poisson initial degree distributions. (b) Derivative of disease prevalence with respect to , i.e., . The jump in with positive values implies a possible discontinuity of . Except for the case , all networks reached their stationary disease prevalence . We here fix the parameters of the system as , and . Each pixel corresponds to the mean computed over 30 Monte Carlo simulations.

Finally, we study the effect of on the observed bifurcations while varying the fraction of initial infections . Instead of using a fixed value of (we recall in Figs. 14–16), we consider to be 0.001, 0.01, 0.2, 0.4, 0.6, 0.8, 0.99, and 0.999. We run simulations to and once again observe stationary states for . In Fig. 17 we show that the AME predicts the disease prevalence curves for the three different initial degree distributions. Figure 17 also illustrates the effect of on the system. For example, in the degree regular case, the infection becomes endemic only for , whereas in the truncated power law case every leads to an endemic state. In contrast, the Poisson case is more involved: The system appears to be bistable, with larger yielding an endemic infection and smaller typically reaching a disease-free state. Even so, we note that the points in the figure include points on both the disease-free and endemic branches for the Poisson case.

FIG. 17.

Bifurcation diagrams of the disease prevalence versus on networks with initial degree distributions that are (a) Poisson, (b) truncated power law, and (c) degree regular. The parameters of the system are , and . Lines are the AME predictions and symbols are the means from 30 simulations. We run simulations for each value , 0.01, 0.2, 0.4, 0.6, 0.8, 0.99, and, 0.999 and plot them with different colors. We set , with all parameters with not close to 1 reaching their stationary states by this time.

We further investigate this behavior with calculations at , so on average we expect to initially see only one infected node, and with probability the networks start with no infected nodes. We performed 30 Monte Carlo simulations for this , finding that 80% of the networks lead to a disease-free state for the Poisson initial condition; on the other hand, 60% of the networks lead to a disease-free state for the truncated power law initial condition. We observe that when is small enough, the system can converge to either an endemic or disease-free state. The AME is unable to predict this bistability from a fixed . However, in the cases where the system exhibits only one type of stationary state, then our AME system gives a good prediction of the stationary disease prevalence . In the degree regular case [see Fig. 17(c)], there is an additional horizontal line indicating that the AME can capture the state when for . In Fig. 17 there is greater variability near because the system does not converge to a stationary state for the simulated times. Figure 17 also shows that for large enough becomes independent of the initial degree distribution, with the results falling on the same universal curve.

DISCUSSION

We have introduced a model variant for the spread of disease on a coevolving network with reinforcement of transitivity, a quintessential property of social networks. Transitivity in the model is controlled by a combination of the rewiring rate and an additional parameter for the probability of closing triangles. This model provides an opportunity to study the role of transitivity in altering the dynamics of disease spread.

We explored the parameter space of and with three different initial degree distributions. We have identified that higher values of lead to lower disease prevalence. In other words, increased reinforcement of transitivity can decrease the disease prevalence. This finding suggests a possible general mechanism for controlling SIS-type epidemics, intuitively encouraging healthy individuals to focus their contacts to healthy individuals within a close social circle, and this will typically do better than building contacts to random people, in the sense of leading to lower endemic rates. The viability of this general intuition beyond the specific setting of direct closure of triangles may be important for future study.

We carried out a bifurcation analysis to understand the properties of the systems at equilibrium. Figure 5 and the Appendix A 2 imply that the phase transition from disease-free to epidemic states as varies might be discontinuous when and continuous when . Validating the nature of this change is a promising potential direction for further studies.

We extended the AME method to include the effect of transitivity reinforcement. We showed that for , our AME system predicts the disease prevalence in the stationary state. Furthermore, we illustrated the accuracies of the stationary degree distributions predicted by the AME. This success of the AME further supports its use to study a variety of binary state dynamics on coevolving networks, accurately predicting properties of large networks at a manageable computational cost. We remark that the high AME accuracy observed here for nonzero may be surprising, in that a key assumption of the AME method is that the networks are locally treelike, but this assumption will seemingly become less valid as local clustering increases with increasing . Nevertheless, our AME method appears able to provide good predictions; the possible reasons for such accuracy remain a potential area for future exploration. (See, e.g., Ref. [41] for other settings where theoretical models remain good predictions in the presence of significant local clustering.)

The behavior of this model at (i.e., only rewiring to close triangles) is different from the cases. When , the system evolves very slowly and does not converge to a stationary state in our simulated times, with the disease prevalence appearing to decay as a small inverse power of . We surmise that the system is slowly approaching a disease-free state. In this case, the AME is not able to quantitatively describe the temporal evolution of the system; however, it appears to still correctly predict the final disease-free state.