key: cord-0004094-t0ma3duq authors: Lee, Hsuan-Wei; Malik, Nishant; Shi, Feng; Mucha, Peter J. title: Social clustering in epidemic spread on coevolving networks date: 2019-06-04 journal: nan DOI: 10.1103/physreve.99.062301 sha: edbcc04de29555ae173e9760dde3575b2718f950 doc_id: 4094 cord_uid: t0ma3duq 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. 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] [6] [7] [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] [11] [12] [13] [14] , and there is extensive literature exploring the impact of network structure on the dynamics of epidemic spread on static networks [10] [11] [12] [13] [14] [15] [16] [17] [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] [23] [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] [26] [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-infectedsusceptible (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 infectedsusceptible 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. We employ the classical setting of the SIS model [30, 31] . We consider a network with N nodes and M undirected edges (both constant), with M = k N/2, where k is the average degree. At a given time, each node is either in a susceptible (S) or infected (I) 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 1 − η) the rewiring is made to another susceptible node selected uniformly at random from the network outside its immediate neighborhood [that is, selfloops 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. Unless stated otherwise, our numerical results of the adaptive SIS system presented here are for networks with N = 25 000 nodes and M = 25 000 edges (that is, mean degree k = 2), 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 1/N, 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 N; a truncated power law distribution with τ = 2.161 and k c = 20 in order to make the mean degree k = 2 (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 k 0 = 2. We observe that at large enough time t, 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 I (expressed as a fraction of nodes) between two consecutive integer times differs by at most 10 −5 . Under this definition, we find that in practice networks with the above three initial degree distributions reach their stationary states before t = 5000 for η < 1. In contrast, when η = 1 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 ). 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 C = (triangles) and the total number of triplets centered at vertex i, respectively [35] . The second measure is the transitivity T (often called global clustering coefficient) [36] , defined as T = 3T /τ , where T 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 C and the transitivity T 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 η = 0, there is no reinforcement of transitivity, i.e., C ≈ 0 (T ≈ 0) 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 η = 1, where we do not observe this convergence, and C (as well as T ) appears to continue to grow with a diminishing rate. We treat η = 1 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 η < 1. Finally, in Fig. 2 we examine the local clustering coefficient C grouped by susceptible nodes S and infected nodes I at time t = 5000 for networks with Poisson degree distributions initially and with different values of η. As η increases, susceptible nodes typically achieve higher values of local clustering coefficient C compared to the infected nodes. Given that our model demonstrates nontrivial transitivity, we now explore how the reinforced transitivity influences the evolution of disease prevalence I. In Fig. 3 we plot the observed disease prevalence at long times in the (γ , η) parameter space, fixing the infection rate β = 0.04 and recovery rate α = 0.005 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 γ = 0, there is no rewiring process, the network does not change, and hence η has no effect. A more interesting behavior occurs for γ 0.05, with large values of disease prevalence that depend FIG. 3 . Phase diagrams in parameter space (γ , η) for the observed disease prevalence I at time t = 10 000 on networks with initial degree distributions that are (a) Poisson, (b) truncated power law, and (c) degree regular. We set here β = 0.04, α = 0.005, and = 0.1. At the selected time t = 10 000 all cases except η = 1 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. on η in a complex manner. Since this parameter regime is most interesting for our present study, in the rest of our simulations we fix γ = 0.04 and vary η to investigate the influence of transitivity reinforcement. 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 (N s ) and the fraction of nodes in the largest component (s 1 ) at t = 5000 for the three initial settings of the degree distribution. We fix the mean degree k = 2, γ = 0.04, β = 0.04, α = 0.005, and = 0.1 and then compute N s and s 1 at t = 5000. Varying η from 0 to 1, we observe the increase in N s and decrease in s 1 , 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 η = 1. 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 η < 1, for all three initial distributions there exists a small bistable region near β . = 0.16. The sudden jump of I ∞ around the critical β when η = 0 (the first column of plots) implies that the transition might be discontinuous when η = 0. This discontinuity appears to shrink as η increases. For η = 0.8, 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 η = 1, with the y axis labeled as I and not I ∞ , highlighting the fact that for η = 1 the system does not reach a stationary state in the observed time. 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 FIG. 5 . Bifurcation diagrams of the disease prevalence in the stationary state I ∞ 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 t = 10 000. The last column belongs to the η = 1 case, noting that the y axis is labeled I (cf. I ∞ ) because the simulations have not reached stationarity. The other parameters of the system are γ = 0.02 and α = 0.005. (We also study the effect of different in initially Poisson networks; see the Appendix A 2.) 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 O(k 2 max ) to O(k 4 max ) equations, where k max 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 S kl (t ) and I kl (t ) be the fraction of susceptible and infected nodes, respectively, of degree k with l infected neighbors at time t [see the diagram in Fig. 6 (a)]. Following the notation in Ref. [34] , we define the zeroth-order moments of the S kl (t ) and I kl (t ) 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, S + I = 1, since the number of nodes is fixed. Similarly, S S + S I + I S + I I = k , because of the conservation of edges. We then have the following ordinary differential equation (ODE) governing the time evolution of the S kl compartment: Similarly, the ODE for the I kl compartment is To describe the AME derivation, we will focus on explaining the S kl equation, as the derivation of corresponding terms in the I kl 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 S dismissing one of its state I neighbors and rewiring to a random node with state S. 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 S and its distance-2 neighbors with state S, breaking links to one of their I 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 S nodes can rewire and be rewired to in this model (whereas I neighbors are only dropped by the action of S nodes). Supposing the center node is of class S (k−1)l , then it has l k−1 proportion of its neighbors infected neighbors, and we describe the number of SI edges among them as (1 − η) , is the effect of a center class with state S rewired to a random node in the network with state S. For the I kl 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 I nodes are rewired. Therefore, the equation lacks terms premultiplied by γ η and γ (1 − η). The resulting system of coupled ordinary differential equations FIG. 6. Illustration of an S node rewiring to a neighbor's neighbor. (a) Before the rewiring, the center nodes in the diagram are of class S kl and I kl , respectively. Consider the two S nodes involved when an SI 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 S (k+1)l . (c) Diagram representing the case when the center actively rewires to a distance-2 neighbor and the center node becomes class S k(l−1) . The other processes of the adaptive SIS system can be similarly diagramed. contains 2(k max + 1) 2 equations, where k max 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 p k (0) is the degree distribution at t = 0. 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. First we plot the disease prevalence I against time t with different η's and different initial degree distributions in Fig. 7 . Note that η = 0 [ 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 I approaches a nonzero value, i.e., the system attains an endemic state. However, in the degree regular case, I 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. 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 t ∈ [0, 1000] for some of the panels of Fig. 8.) 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 η < 0.8. For 0.8 η 1, the AME underestimates the stationary disease prevalence I but still captures qualitative features of the evolution. Remarking again that η = 1 is a special case, we note that the system does not reach a well-defined stationary state in our simulations, with I continuing to decrease. In Fig. 9 we further explore the temporal evolution of I when η = 1 while varying k from 2 to 10. We identify a slow power law decay in I for k = 2; using the curve fitting in MATLAB ® , we find that the disease prevalence I for η = 1 and k = 2 scales as I = 1.9t −0.173 (also shown in Fig. 9 ), with I appearing to slowly approach a disease-free state. As the mean degree increases, the observed decay with t 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 η = 1 is qualitatively very different, it does appear to predict a correct I ∞ = 0 disease-free stationary state (see Fig. 8 ). To summarize, we plot simulation results (markers) and AME-based estimates (lines) of I ∞ 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.) 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 η = 0.4 for the three different initial distributions; that is, the degree distributions of the susceptible S and infected I 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 SI edges in the system; however, there are still a few II edges and isolated infected nodes and as a result we still see infected nodes with different degrees in the stationary states. To determine how η impacts the degree distribution, we plot the initial and stationary degree distributions for different η's in Fig. 12 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 k = 2 [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. Similar to Fig. 5 in Ref. [34] , we plot (Fig. 13 ) bifurcation diagrams of stationary disease prevalence level I ∞ versus infection rate β with an initial Poisson degree distribution and with different mean degrees k = 2, 4, and 8. The other parameters here include γ = 0.02 and α = 0.005. We obtained the average stationary disease prevalence level I ∞ at each value of β from 30 Monte Carlo simulations for each initial infected proportion = 0.001, 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 k = 2, from 0 to 0.02 in Figs. 13(e)-13(h) for k = 4, and from 0 to 0.008 in Figs. 13(i)-13(l) for k = 8. As evident in the figure, our AME calculation reasonably predicts the stationary disease prevalence level I ∞ , including its phase transition and the bistability near that transition. 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 η = 1 (see Fig. 9 ). To further investigate whether the nature of these transitions might be a finite-size effect, in Fig. 14 2 )/ I ∞ , as defined in [40] , for different network sizes N. With the peak of the susceptibility identifying the transition, the results in the figure show that for η < 1 the network size N does not greatly change the transition location in β, except in the smallest networks. However, we again see a qualitative difference for η = 1, making a definitive conclusion in this case more difficult. 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 N = 10 000 to N = 50 000 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 k = 4 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 I ∞ 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 N, β, 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. Furthermore, to demonstrate the influence of η on the transitions, in Fig. 16(a) we plot a phase diagram for the observed stationary disease prevalence I ∞ 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, dI ∞ /dβ, observing only a narrow range with a large positive value near the transition points apparent in Fig. 16(a) . Also, the jump in dI ∞ dβ occurs at a lower β as η increases, again showing the dependence of the transitions on η. Finally, when η = 1, we see no jumps of I ∞ when β varies, further suggesting that the system behavior at η = 1 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. 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 = 0.1 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 t = 10 000 and once again observe stationary states for η < 1. 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 > 0.1, 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 diseasefree state. Even so, we note that the = 0.001 points in the figure include points on both the disease-free and endemic branches for the Poisson case. We further investigate this behavior with calculations at = 0.000 04, so on average we expect to initially see only one infected node, and with probability (1 − 0.000 04) 25 000 . = 37% 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 I ∞ . In the degree regular case [see We vary η from 0 to 1, recalling that η = 1 does not reach stationarity by t = 20 000. We note that I ∞ is consistent across network sizes in the region tested in this paper. Fig. 17(c) ], there is an additional horizontal line indicating that the AME can capture the state when I ∞ = 0 for 0.1. In Fig. 17 there is greater variability near η = 1 because the system does not converge to a stationary state for the simulated times. Figure 17 also shows that for large enough , I ∞ becomes independent of the initial degree distribution, with the results falling on the same universal curve. 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 diseasefree to epidemic states as β varies might be discontinuous when η < 1 and continuous when η = 1. 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 η < 0.8, 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 η = 1 (i.e., only rewiring to close triangles) is different from the η < 1 cases. When η = 1, 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 t. We surmise that the system is slowly approaching a diseasefree 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. Different cutoff degrees in the truncated power law distribution have different levels of disease prevalence temporarily; however, as t increases, they lead to the same level of disease prevalence. Here we revisit the bifurcation diagram in Fig. 5 and let the initial infected ratio vary from 0.001 to 0.999. In Fig. 19 we can see that the critical region of transition also depends on . To investigate how the AME describes the temporal evolution of the system dynamics, especially the transitory dynamics in the early periods, we first present in Fig. 20 a zoomed-in version of Fig. 7 for the region t ∈ [0, 1000]. Overall, the AME is able to reproduce the general shape of the time trajectories of the transitory dynamics. In particular, we note that the AME performs well when FIG. 19 . Bifurcation diagrams of the disease prevalence in the stationary state I ∞ versus β at various values of η and . Each column of plots corresponds to a particular value of η and each row of plots corresponds to a particular value of . As in Fig. 5 , every point represents a mean of 30 realizations of the system at t = 10 000. The last column belongs to the η = 1 case, noting that the y axis is labeled I (cf. I ∞ ) because the simulations have not reached stationarity. The networks have a Poisson degree distribution initially. The other parameters of the system are γ = 0.02 and α = 0.005. initially the network is degree regular, regardless of η. For the other two initial degree distributions, the AME performs well when η is not too large, with the discrepancy between simulations and AME approximations becoming larger as η increases, especially for the truncated power law cases. Moreover, inspired by Ref. [34] , in Fig. 21 we include the evolution of the fraction of SI links S I , the effective branching factor κ S IS ≡ S SI /S I , and the average number of connections that susceptible nodes share with other susceptibles C SS ≡ S S /(S S + S I ) as well as the AME approximation for these quantities. In an SIS model, S I measures the level of links that could potentially pass the disease, κ S IS measures the average number of susceptible neighbors that the infected end of a SI link has, and C SS measures the fraction of susceptible neighbors of a susceptible node. In Fig. 21 we observe that the AME captures the qualitative behavior of the temporal evolution of the different quantities, but again the discrepancy between simulations and AME approximations becomes larger as η increases. Having shown at multiple places in the main text that the AME accurately predicts the stationary level of disease prevalence and various quantities of interest, we also note here the level of agreement between the AME and simulation for the stationary levels of the quantities in Fig. 21 , with good agreement at η = 0 and increased error for η = 0.4. On the other hand, the differences between simulations at different η are not particularly large for η up to around 0.6 and with qualitative differences only for η = 1 (see Figs. 8 and 10). The changes in the stationary levels predicted by the AME do not vary greatly with η either, but they appear to do so in the correct direction for these quantities here. While the AME does not capture the system dynamics well when η is large, it nevertheless provides good predictions for the stationary disease prevalence (Figs. 10, 13, and 17) and degree distributions (Figs. 11 and 12) for η not too close to 1. Networks and epidemic models Spread of epidemic disease on networks Adaptive coevolutionary networks: A review Modeling complex systems with adaptive networks Impact of Non-Poissonian Activity Patterns on Spreading Processes Small but slow world: How network topology and burstiness slow down spreading How memory generates heterogeneous dynamics in temporal networks Birth and death of links control disease spreading in empirical contact networks Behavioral and neural correlates to multisensory detection of sick humans Dynamic social networks and the implications for the spread of infectious disease Infectious disease modeling of social contagion in networks Fluctuating epidemics on adaptive networks Infection spreading in a population with evolving contacts Modelling the influence of human behavior on the spread of infectious diseases: A review Epidemics and percolation in small-world networks Percolation and epidemics in random clustered networks Small World Effect in an Epidemiological Model Forecast and control of epidemics in a globalized world Random Graphs with Clustering Bond percolation on a class of clustered random networks Epidemic spreading on evolving signed networks Dynamics of social networks Nonequilibrium phase transition in the coevolution of networks and opinions The effect of a prudent adaptive behavior on disease transmission Role of social environment and social clustering in spread of opinions in coevolving networks Transitivity reinforcement in the coevolving voter model Dynamical origins of the community structure of an online multi-layer society High-Accuracy Approximation of Binary-State Dynamics on Networks Binary-State Dynamics on Complex Networks: Pair Approximation and Beyond Mathematical Models in Population Biology and Epidemiology Epidemic Dynamics on an Adaptive Network Largenet2: An object-oriented programming library for simulating large adaptive networks A general method for numerically simulating the stochastic time evolution of coupled chemical reactions Adaptive networks: Coevolution of disease and topology Collective dynamics of "smallworld" networks Social Network Analysis: Methods and Applications Graph fission in an evolving voter model Link-based formalism for time evolution of adaptive networks Evolutionary prisoner's dilemma games coevolving on adaptive networks Epidemic thresholds of the susceptible-infected-susceptible model on networks: A comparison of numerical and theoretical results The unreasonable effectiveness of tree-based theory for networks with clustering The research reported herein was supported by the Eunice Kennedy Shriver National Institute of Child Health and Human Development of the National Institutes of Health under Award No. R01HD075712. Additional support was provided from the James S. McDonnell Foundation 21st Century Science Initiative, Complex Systems Scholar Award Grant No. 220020315. The content is solely the responsibility of the authors and does not necessarily represent the official views of the institutions supporting this work. Truncated power laws depend on an additional parameter, known as the cutoff degree, k c [see Eq. (2) above for details].Here we examine the effects of k c on the variable I (disease prevalence level) at time t (see Fig. 18 ). Specifically, we keep the mean degree k = 2 and vary the cutoff degree k c from 20 to 60. We found that a different cutoff degree of the truncated power law distribution yields a different level of disease prevalence level temporarily, but they have the same level of disease prevalence and degree distribution at stationarity.