# Synergistic interactions promote behavior spreading and alter phase transition on multiplex networks

###### Abstract

Synergistic interactions are ubiquitous in the real world. Recent studies have revealed that, for a single-layer network, synergy can enhance spreading and even induce an explosive contagion. There is at the present a growing interest in behavior spreading dynamics on multiplex networks. What is the role of synergistic interactions in behavior spreading in such networked systems? To address this question, we articulate a synergistic behavior spreading model on a double layer network, where the key manifestation of the synergistic interactions is that the adoption of one behavior by a node in one layer enhances its probability of adopting the behavior in the other layer. A general result is that synergistic interactions can greatly enhance the spreading of the behaviors in both layers. A remarkable phenomenon is that the interactions can alter the nature of the phase transition associated with behavior adoption or spreading dynamics. In particular, depending on the transmission rate of one behavior in a network layer, synergistic interactions can lead to a discontinuous (first-order) or a continuous (second-order) transition in the adoption scope of the other behavior with respect to its transmission rate. A surprising two-stage spreading process can arise: due to synergy, nodes having adopted one behavior in one layer adopt the other behavior in the other layer and then prompt the remaining nodes in this layer to quickly adopt the behavior. Analytically, we develop an edge-based compartmental theory and perform a bifurcation analysis to fully understand, in the weak synergistic interaction regime where the dynamical correlation between the network layers is negligible, the role of the interactions in promoting the social behavioral spreading dynamics in the whole system.

## I Introduction

A central problem in network science and engineering is to understand, predict, and control the dynamics of virus or information spreading on complex networks Barrat et al. (2008); Vespignani (2012); Pastor-Satorras et al. (2015). Social contagion processes such as the propagation of an opinion, diffusion of a belief, and spread of a particular behavior, occur commonly in the real world Valente (1996); Christakis and Fowler (2007); Young (2009); Rogers (2010); Centola (2011); Banerjee et al. (2013); Zha et al. (2016); Zhang et al. (2016). With the modern technological advances, a variety of online social networking platforms (e.g., Facebook and Youtube) have become a routine necessity for a substantial fraction of individuals in the entire population. Spreading dynamics in modern online social networks have attracted a great deal of recent attention and a variety of mathematical models have been articulated to understand and predict the relevant phenomena Pastor-Satorras et al. (2015); Castellano et al. (2009); Dodds and Watts (2004a, 2005). For example, the threshold model, a binary state spreading model, was introduced earlier to address the phenomenon of behavior adoption, where a node in a social network adopts a new behavior only when the number Granovetter (1973) or the fraction Watts (2002) of its nearest adopted neighbors exceeds a threshold value. A representative threshold model reveals the phenomenon that the final size of the nodes adopting the behavior first grows continuously and then decreases discontinuously as the mean degree of the network is increased Watts (2002). Within the threshold model, the effects of parameters and network structure on the dynamics of social behavioral spreading have been studied, which include the initial seed size Gleeson and Cahalane (2007), the clustering coefficient Whitney (2010); Hackett et al. (2011); Zhuang et al. (2017), the community structure Gleeson (2008); Nematzadeh et al. (2014) and multiplexity Brummitt et al. (2012); Yağan and Gligor (2012); Cozzo et al. (2013). The dynamical process described by the threshold model, however, is Markovian because the state of a node depends only on the current state of its neighbors. The original model is thus not able to encompass an important aspect of real contagion dynamics: social reinforcement originated from the memory effect Dodds and Watts (2004b); Krapivsky et al. (2011); Zheng et al. (2013); Liu et al. (2016) - a feature that is characteristically non-Markovian. To overcome this deficiency of the classical threshold model, a non-Markovian behavior spreading model taking into account the received cumulative pieces of behavioral information for any node to adopt the behavior was introduced Wang et al. (2015a). A prediction of the modified model is that the dependence of the final behavior adoption size on the information transmission rate can change from being discontinuous to being continuous through continuous changes in the dynamical or structural parameters. The non-Markovian behavior spreading model also allows additional issues such as the heterogeneity of adoption thresholds Wang et al. (2016), the limited contact capacity Wang et al. (2015b), and the effect of temporal network structure Liu et al. (2017a) to be addressed.

Most previous works on network behavior spreading focused on a single social behavior contagion process through empirical methods Centola (2011); Banerjee et al. (2013) and mathematical models Dodds and Watts (2004a, 2005); Granovetter (1973); Watts (2002); Wang et al. (2015a); Liu et al. (2017b). In the real world, it is common for two or more distinct behaviors to spread simultaneously in a social system, where interactions between the corresponding spreading processes inevitably arise. For example, individuals who have adopted Windows services are more likely to use other services from the same company, e.g. Microsoft Office. In online networking systems, two different tweets on the same event or subject can diffuse on the twitter network at the same time. The user seeing one tweet will experience an increased exposure to the other tweet, and vice versa, since these two tweets are closely related. In this case, the two tweets spread synergistically as they mutually prompt each other in the process of retweeting Myers and Leskovec (2012). The synergistic mechanism is also typical in the adoption of online services. A good example is the adoption of two online services, say Google and Youtube through two types of tweets: one containing the URLs with google and another with youtube. The numbers of the two types of tweets are synchronized most of the time, implying that they are synergistic to each other Zarezade et al. (2015). The synergistic effect also occurs in disease spreading, where the interaction between pathogens may mutually strengthen their spreading process, and such an effect may have played a role in the co-epidemic of the Spanish flu and pneumonia in 1918 Taubenberger and Morens (2006); Brundage and Shanks (2008); Chen et al. (2013); Cai et al. (2015); Hébert-Dufresne and Althouse (2015). In spite of its ubiquity, the synergistic mechanism among two or more simultaneously spreading behaviors was not investigated in previous studies Dodds and Watts (2004a, 2005); Granovetter (1973); Watts (2002); Wang et al. (2015a).

In this paper, we articulate a synergistic social behaviors spreading model to address and understand the impacts of synergistic interactions among multiple behaviors on their spreading. As the spreading of each behavior typically occurs on a different network layer, it is necessary to incorporate a multilayer network structure Domenico et al. (2016); Kivelä et al. (2014); Boccaletti et al. (2014). To be concrete, we consider the spreading dynamics of two distinct behaviors in two-layer coupled networks, where each layer supports the spreading of one behavior with its own transmission path, as described by a non-Markovian process. The synergistic mechanism between the two behavior adoption dynamics is that, once a node adopts a behavior in one layer, it becomes more susceptible to adopting the other behavior that spreads in the other network layer. We develop an edge-based compartmental theory to analyze and understand how the synergistic interactions impact the simultaneous spreading dynamics of the behaviors. We find, as suggested by intuition, that the synergistic interactions greatly facilitate the adoption of both behaviors. However, surprisingly, a phenomenon is that the adoption of one behavior can lead to a characteristic change in the adoption of the other behavior: its final adoption size versus its information rate can change from being discontinuous to continuous, where the former corresponds to a first-order phase transition while the latter to a second-order transition. Remarkably, the synergistic effect can induce a two-stage contagion process, in which nodes having adopted one behavior in one layer will adopt the other behavior in the other layer. When there is a sufficient number of seeds, i.e., when the number of nodes having adopted the other behavior in the other layer is sufficiently large, the remaining nodes will adopt the behavior quickly. While it is intuitively understandable that the synergistic interactions can promote the spreading dynamics of the distinct behaviors involved, our work lays a quantitative foundation for this phenomenon. Our model will not only serve as a useful framework to understand the interplay between synergy and simultaneous spreading of multiple behaviors or diseases, but will also provide insights into predicting or even controlling the underlying dynamics. Due to the ubiquity of synergy in different fields such as social science, computer science, biology and biomedicine, broad relevance of our model is warranted.

In Sec. II, we describe the network and the synergistic behavior spreading models. In Sec. III, we carry out a detailed theoretical analysis. In Sec. IV, we present extensive simulation results with respect to the theoretical predictions. In Sec. V, we summarize the main results and discuss a few pertinent issues.

## Ii Model

There are two components in our model: multiplex networks and spreading dynamics of synergistic behaviors. We first introduce the model of multiplex networks, and then present the synergistic behavior spreading model.

### ii.1 Model of multiplex networks

In general, network layers in an interdependent networked system have different internal structures and dynamical functions. To capture the essential dynamics of simultaneous spreading of distinct behaviors, we focus on multiplex networks Domenico et al. (2016); Kivelä et al. (2014); Boccaletti et al. (2014). Consider the simple setting of a duplex system consisting of two layers or subnetworks. Initially, we generate two independent layers, denoted as and , which have the same node set and support the spread of behaviors and , respectively. We use the configuration model Catanzaro et al. (2005) to generate each subnetwork, where the degree distribution of layer is completely independent of the distribution of layer . For large and sparse subnetworks, the configuration model stipulates that both interlayer and intralayer degree-degree correlations are negligible.

### ii.2 Synergistic behavior spreading model

We use a representative non-Markovian spreading model, the susceptible-adopted-recovered (SAR) Wang et al. (2015a) model, to describe the dynamics of behavior spreading, and then introduce the synergistic mechanism between the spreading processes of the two behaviors.

For each behavior , at any time a node will be in one of the three states: susceptible (), adopted () and recovered (). A node in state has not adopted behavior but it has an interest in . A node in the state has adopted the behavior and can transmit the information about the behavior (denoted as information ) to its neighbors. The node loses interest in transmitting the information when it is in the state. The evolution process of behavior can be described, as follows. Initially, fraction of nodes are randomly chosen as the nodes that have adopted the behavior and the remaining nodes are set to be in the susceptible state. At each time step, each node in the state transmits the information to each of its susceptible neighbors with the transmission rate . Suppose a neighboring node already has accumulated pieces of information from its distinct neighbors. One more successful transmission will make the number of information pieces to become . We assume non-redundant information transmission, i.e., once an adopted node has transmitted the information to node , the former will not transmit the same information to latter again. If the cumulative number pieces of information that the susceptible node has is equal to or larger than a threshold, the node will adopt the behavior and changes its state to . Simultaneously, each node will turn to the state at the recovery rate . The behavior spreading process will terminate when all the adopted nodes have recovered. More specifically, and are the fractions of nodes randomly chosen as seeds (i.e., adopted nodes) for behavior and on each layer, respectively, where the remaining nodes are in the susceptible state. Information () diffuses in layer () with transmission rate (), and the recovery rates for behaviors and are and , respectively.

In the general SAR model, each susceptible node has its own adoption threshold for a behavior. However, for simplicity in modeling the synergistic interaction between the spreading of the two behaviors, we assume that all nodes have the same adoption threshold for each behavior: we denote the adoption threshold for behavior 1 in layer as and that for behavior 2 in layer as . As a manifestation of mutual synergy, a node having adopted one behavior will become more susceptible to adopting the other behavior. To quantify the synergistic effect, we assume that, once node has adopted behavior (), it will generate an increase () in the number of pieces of information about behavior (). The quantities and thus characterize the strength of the synergistic effect, and we have and . For , a node having adopted behavior in layer will not impact on its adoption of behavior in layer . Similarly, the adoption of behavior will have no effect on adopting behavior if . If a node has adopted behavior , it will adopt behavior only if , where represents the number of cumulative pieces of behavioral information in layer that this node has received from distinct neighbors.

## Iii Theory

We exploit the edge-based compartmental theory Wang et al. (2015a); Miller et al. (2012); Yang and Zhou (2012); Karrer and Newman (2010) to analyze the dynamical process of behavior spreading subject to synergistic interactions, under the assumption that each subnetwork is large and sparse with no internal degree-degree correlations. We also assume that the degree distribution of network is completely independent of that of network , so interlayer degree-degree correlation can be neglected too. The fraction of nodes in each state can be treated as a continuous variable. For each behavior , we denote , and as the fractions of nodes being in the susceptible, adopted, and recovered state, respectively, for behavior in the corresponding layer at time . During the spreading process, the susceptible nodes adopting behavior decreases the value of but leads to an increase in , and the recovery of the adopted nodes for behavior decreases but increases . Using these notations, the dynamical evolution equations for behavior can be written as

(1) |

and

(2) |

For , the states of all individuals remain unchanged and is the final adoption fraction of behavior .

### iii.1 Edge-based compartmental theory

Despite that the spreading processes of behaviors 1 and 2 occur in different networks ( and , respectively) and the dynamical parameters such as the information transmission rates ( and ), the recovery rates ( and ), and the adoption thresholds ( and ), are different, the mathematical equations governing the underlying processes have identical forms. It thus suffices to derive the equations for behavior spreading in layer .

To solve Eqs. (1) and (2), we need to calculate the fraction of susceptible nodes for behavior at time step . Firstly, for nodes of degree in layer , two cases can arise where the nodes do not adopt behavior 1: (1) these nodes have not adopted behavior on layer and the cumulative number of received pieces of information in layer is less than , and (2) these nodes have already adopted behavior in layer , but the cumulative number of received pieces of information in layer is less than . Under the assumption that there is no dynamical correlation between the layers, we have that the fraction of susceptible nodes of degree for behavior at time is given by

(3) | |||||

In Eq. (3), the first term on the right side is the probability that a node of degree in layer at time does not adopt behavior 1. This term contains two parts that describe the following two situations, respectively: (1) the received cumulative number of pieces of information 1 is less than with probability , and (2) with probability , a random node in layer does not adopt behavior 2 at time (i.e., a node in layer does not adopt behavior 2 and is still in the susceptible state), where the quantity is the probability for a node of degree to have received pieces of information by time in layer . Combining the two parts, we find that the first term is identical to the second term in Eq. (3). Using the degree distribution of network , we can express the fraction of susceptible nodes for behavior as

(4) |

In Eq. (3), the quantity can be expressed as

(5) |

where denotes the binomial distribution and is the probability that a random neighbor of node in layer has not transmitted the behavioral information to node by time . To take into account the dynamical correlations among the states of the adjacent nodes, we make use of the cavity theory Wang et al. (2015a); Miller et al. (2012); Yang and Zhou (2012); Karrer and Newman (2010) to analyze the quantity , where node is in the cavity state so that it cannot transmit the behavioral information to its neighbors but it can receive the information from its neighbors.

To solve Eqs. (3) and (4), we need the value of [the computation of is the same as that of ]. Noting that a random neighbor of node in layer can be in one of the following three states: , and , we have that is the sum of the probabilities that the neighbor does not transmit information to when is in the , or state. We have

(6) |

where [ or ] denotes the susceptible (adopted or recovered) neighbor of which has not transmitted information to node up to time in layer .

Suppose a random neighbor of degree of node is susceptible initially, node cannot transmit information to since is in the cavity state. Node can only receive the information from its other neighbors. The probability that node has received pieces of information in layer by time is then

(7) |

Similar to Eq. (3), we have that the probability that the neighboring node is still in the susceptible state for behavior at time is given by

For uncorrelated networks, the probability for a random edge to connect a node of degree is , where is the average degree of network layer . A neighboring node in the susceptible state cannot transmit the behavioral information. Thus, is equal to the probability that the neighboring node is in the susceptible state, which is

(9) |

If a random neighbor is in the adopted state for behavior 1, success in information transmission from node to node will result in a decrease in . We thus have

(10) |

At the same time, once the adopted neighbor has recovered before it can transmit information to node , there will be an increase in . (Note that here we use the synchronous updating rule, meaning that the transmission and recovery events happen consecutively in discrete time steps.) The increase in contains two parts that describe the following two situations, respectively: (1) with probability , the neighboring node has not transmitted information 1 to , and (2) simultaneously, node recovers with probability . Combining these two parts, we obtain the increment of as

(11) |

Combining Eqs. (10) and (11), we obtain an explicit expression for :

(12) |

Inserting Eqs. (9) and (12) into Eq. (6), we can write as

(13) | |||||

Substituting Eq. (13) into Eq. (10), we get the time evolution of as

(14) | |||||

Following a similar procedure, we can derive the expression of , the probability that a random neighbor of node in layer has not transmitted the behavioral information to node by time , and . We have

(15) | |||||

and

(16) | |||||

where the form of in Eq. (15) is similar to , and in Eq. (16) is similar to . It is thus not necessary to write down the expressions again. Using the degree distribution of network , we have the fraction of susceptible nodes at time in layer as

(17) |

Iterating Eqs. (1)-(4) and (14)-(17), we can obtain the fractions of susceptible nodes at time in both layers: and . In addition, we can substitute [] into Eqs. (1) and (2) and calculate the fractions of the adopted nodes and of the recovered nodes in layer () at time . Taking the limit , we can obtain the final fractions of adoption of the two behaviors. Results on the final adoption fractions from direct numerical simulations together with the corresponding theoretical predictions for different parameter values are shown in Fig. 1. We obtain a good agreement between theory and numerics. For example, for and , Fig. 1(b) shows that, without the synergistic effect of behavior 1, i.e., , behavior 2 will not exhibit any outbreak. For , behavior 2 is adopted globally. When there are mutual synergistic effects, e.g., and or and , the adoption of both behaviors is enhanced, as shown in Figs. 1(c) and 1(d), respectively. Note that there are some outliers (e.g., there are one black square in Fig. 1 (a) and two black squares in Fig. 1 (d)) around the critical transmission rate since the SAR model is not a deterministic threshold model, which is in contrast to the Watts threshold model. The randomness exists in the process of simulations when the behavior information transmission rate is smaller than 1. Supposing a susceptible node with adoption threshold equal to 3, when it has three adopted neighbors it will not adopt the behavior if one of its adopted neighbor does not succeed in transmitting the behavior information. As shown in the inset of Fig. 1(d), there are some stochastic simulations that does not increase from a very smaller value to a value close 1 directly.

A fundamental issue in spreading dynamics in complex networks is phase transitions Pastor-Satorras et al. (2015). As a system parameter (e.g., the infection rate) changes through a critical point, the final size of the infected nodes starts to increase from zero. An abrupt and discontinuous increase in the final size signifies a first-order phase transition, while a gradual and continuous change is indicative of a second-order phase transition. An objective of our study is then to uncover and understand the effect of synergistic interactions on the phase transitions associated with the social behavior spreading dynamics. To analyze the phase transition, we focus on the fixed point (root) of Eqs. (14) and (15) associated with the final state (i.e., ). Simplifying notation as and , we write Eqs. (14) and (15) as

(18) |

and

(19) |

respectively, where

(20) | |||||

and

(21) | |||||

Because of the nonlinear functions in Eq. (20) and in Eq. (21), to analyze the whole parameter space is infeasible. We thus focus on some representative or benchmark cases to gain certain analytic understanding of the numerical results. Specifically, we consider two cases in terms of the adoption thresholds of the two behaviors: (1) the adoption threshold of one behavior is less than that of the other behavior ( or ), and (2) .

### iii.2 Solutions for

For , and , indicating that the adoption of behavior has no effect on the spread of behavior but the adoption of the latter will enhance the spread of former, as shown in Fig. 1. Because Eqs. (18) and (19) are nonlinear functions of and , typically there are multiple roots. In addition, there is persistent transmission of behavioral information from individuals in an adopted state (i.e., or ) to their neighbors, so and decrease with time. Thus, if Eqs. (18) and (19) possess more than one stable fixed point, only the one with the maximum value is physically meaningful Wang et al. (2015a). Since Eq. (18) contains the parameters and only, for a given value of , we can obtain the value of . For given values of the parameters and , with we can solve Eq. (19) numerically. As shown in top panel of Fig. 2, we see that Eq. (19) typically has a non-zero trivial solution even for small values of , indicating that, even when the initial adopted fraction of behavior is small (e.g., ), it will always be adopted by a certain fraction of the nodes. However, the initial fraction of seeds will have an effect on the final adoption size Gleeson and Cahalane (2007); Wang et al. (2015a). To better focus on the effect of synergistic interactions on simultaneous spreading of the two behaviors, we set and calculate the final adoption size versus the behavioral information transmission rate with a particular eye on the possible type of phase transitions.

For , the number of roots (fixed points) of the function is 1 or 3, as shown in Fig. 2(a). Because the physically meaningful solution is the maximum value of the stable fixed point of Eq. (19), there is no global outbreak in behavior [verified numerically, see Fig. 4(a)]. For , the function is tangent to the horizontal axis at for the critical value of . Further increasing above removes the tangent point and leaves with only one intersection point with the horizontal axis. Importantly, from the standpoint of bifurcation analysis, we see that, at this point, the physically meaningful fixed point decreases abruptly to a small value, signifying a first-order phase transition. The critical value for a given can be obtained by using the criterion that a nontrivial solution of Eq. (19) emerges, which corresponds to the point at which the function is tangent to horizontal axis at the critical value of . That is, the critical condition for this case can be obtained by combining Eqs. (18) and (19) and the following equation

(22) |

For and , Eq. (19) has a single root whose value decreases with , as shown in Fig. 2(c). This means that increases with continuously.

For a given value of the transmission rate of behavior , the critical condition is then that behavior will be adopted if its transmission rate is larger than . Similarly, we can compute the minimal information transmission rate of behavior 1 required for a global outbreak of behavior . In particular, setting to be the maximum value (i.e., ) and substituting it into Eqs. (19) and (22), we get the critical values of and . Substitute these values into Eq. (18), we obtain , the minimal information transmission rate of behavior .

Numerical solutions of Eq. (19) also show that, for large values of and , it has one fixed point only when varying , so increases with continuously. As a result, there exists the critical parameter value (i.e., ), across which the dependence of on changes from being discontinuous to continuous. For the special case of (e.g., , , and ), we can numerically solve Eqs. (19) and (22), together with the condition Baxter et al. (2010)

(23) |

Once is determined, we can substitute the value of into Eq. (18) to get . In particular, increases with discontinuously for and the increasing pattern becomes continuous for . Using the same approach, we can determine the critical value of above (below) which increases with discontinuously (continuously).

### iii.3 Solutions for

This is the symmetric case where , , , and . The symmetry implies and . For simplicity, we denote and . Equations (18)-(21) can be written as

(24) |

where

Similar to treating Eq. (III.1), we have

The final fraction of the susceptible nodes of behavior () in layer () is given by

(26) | |||||

Using the same analysis method as for the case , we find that the number of fixed points of Eq. (24) is or , as shown in the lower panel of Fig. 2. Whether there is a tangent point between the function and the horizon axis depends on the strength of synergistic interactions. For , there is no tangent point and only the maximum value of the fixed point of Eq. (24) is physically meaningful, indicating that behavior is adopted by a small fraction of nodes only. For and , the function can be tangent to the horizon axis, as shown in Figs. 2(e) and 2(f). When is increased passing through , the tangent point disappears and the function has only one intersecting point with the horizontal axis. In this case, the fixed point changes discontinuously to a small value, signifying a first-order phase transition.

## Iv Numerical validation

In this section, we perform extensive simulations of behavior spreading on different multiplex networks. We use the notation “RR-RR” to denote the case where both layer and layer host the random regular networks. The notation “ER-SF” represents the setting where layer is an Erdös-Rényi (ER) random network Erdös and Rényi (1959) and layer hosts an scale-free (SF) network Barabási and Albert (1999). Other possible combinations are “ER-ER”, “SF-SF” and “SF-ER”. The size of each network is and the average degree is for both networks. The initial adoption fractions of behavior in layer and behavior in layer are set to be . To calculate the pertinent statistical averages, we use multiplex network realizations and at least independent dynamical realizations for each parameter setting. Unless otherwise specified, the above parameters are adopted in the simulations. Let denote the situation where a node is in the or state in layer so, for example, the notion means that, in layer , a node is in the adopted state or recovered state but it is in the susceptible state in layer . Similarly, indicates that a node is in the adopted state in layer and is in the susceptible state in layer , which means that the node adopts behavior but not behavior 2.

### iv.1 RR-RR multiplex networks

We first perform direct numerical simulations of behavioral spreading dynamics on double layer networked systems consisting of two random regular networks to provide support for our theoretical predictions.

Our theoretical analysis in Sec. III.2 gives that, for , synergistic interactions can promote behavior adoption and spreading. To be concrete, we set and . Figure 3(a) shows the time evolution of the fraction of the recovered nodes in layer for different values of the synergistic interaction strength . We see that behavior will not outbreak if . For and , exhibits a two-stage contagion process, where nodes having adopted behavior in layer will first adopt behavior , until when there is a sufficient number of seeds (i.e., nodes having adopted behavior ) in layer to stimulate the remaining nodes. When this happens, behavior will be adopted quickly in layer . This phenomenon can be explained by noting that, for a small fraction of the initial seeds for behavior 2 [i.e., ], if the synergistic effect of adoption of behavior 1 is absent [i.e., ], behavior 2 will not be adopted globally and only the recovery of the seeds can lead to an increase in the value of . Note that the number of nodes increases with the adoption of behavior in layer [Fig. 3(b)] since the nodes will change to nodes and there is no decrease in the number of nodes in the network. For , nodes that have adopted behavior 1 are more likely to adopt behavior 2 as compared with those that have not adopted behavior 1. Nodes having adopted behavior 1 in layer will first adopt behavior 2 in layer , as indicated by the decrease in the number of the nodes in Fig. 3(c). Before most of the nodes have adopted behavior 2, the seeds (i.e., adopted nodes for behavior 2) in layer are sufficient to stimulate the remaining nodes to adopt behavior 2, inducing a two-stage contagion process. A similar phenomenon occurs for . When the simulation results are compared with the theoretical predictions, we find the former matches well with the latter for . While the deviation emerges when , which are derived from the finite-size effects of the networks and the dynamical correlation between layers. From the bottom panels of Fig. 3, we will find the deviation is decreased when increasing the network size, but the deviation will still exist since the interlayer dynamical correlations are ignored in the theoretical method.

Figure 4(a) shows, for , and , versus for different values of , where the fraction of the nodes in the system is about . As the synergistic interaction strength is increased, behavior is adopted more readily since the number of information pieces about it is decreased. A remarkable phenomenon is the characteristic change in the dependence of on . In particular, for , increases with discontinuously but the increasing pattern becomes continuous for . The reason for the characteristic change is that, for , the nodes having adopted behavior still need to receive additional two (i.e., ) pieces of information to adopt behavior . The system will accumulate a relatively large number of nodes in the subcritical state when the behavioral information transmission rate approaches the critical point, as shown in the inset of Fig. 4(a). Therein, the subcritical state is defined as the node in such state will adopt the behavior if it receives one additional piece of behavior information Wang et al. (2015a). A slight increase in will cause a node in this state to receive an additional piece of information and thus adopts behavior . The node can then transmit the information to its neighbors, which will cause its subcritical neighbors to adopt behavior accordingly, and so on, leading to an avalanche of behavior adoption for the nodes. When most of the nodes have adopted behavior 2 in an abrupt fashion, there is a sufficient number of nodes in layer to stimulate the remaining nodes to adopt behavior 2. As a result, increasing slightly can lead to a discontinuous change in the value of . However, for , only one additional piece of information about behavior 2 is needed for the nodes to adopt this behavior. As the value of is increased from zero, some nodes may receive one piece of information about behavior 2 and adopt it, leading to a continuous decrease in the number of nodes in the subcritical state, as shown in the inset of Fig. 4 (b). This is equivalent to the dynamical process in the susceptible-infected-recovered (SIR) model, in contrast to the cascading process in, for example, the Watts threshold model. As a result, the value of first increases with continuously. When most of nodes have adopted behavior 2, the fraction of adopted nodes in layer is sufficient to stimulate the remaining nodes to adopt behavior 2. Since the fraction of adopted nodes is relatively large [e.g., ], the value of increases with continuously Wang et al. (2015a) at a faster rate, as shown in Fig. 4(a). The same process occurs for . These numerical results agree well with our bifurcation analysis based theoretical prediction.

Figure 4(b) shows the dependence of on for different values of . For a relatively small value of (e.g., ), increases with continuously, which can be understood by noting that, in this case, a global adoption of behavior 2 requires more seeds in layer , and the spread of this behavior depends strongly on the spread of behavior 1. However, for relatively large values of (e.g., and ), versus can exhibit an abrupt or discontinuous increase. In this case, a slight increase in the fraction of seeds for behavior 2 is sufficient for it to spread globally by its own dynamics. Both the continuous growth for small values of and the discontinuous increase for larger values of are predicted by our bifurcation analysis based on Eqs. (18), (19), (22) and (23) by replacing with in Eqs. (22) and (23). There is a good agreement between numerics and theory.

Our analysis and numerical computations indicate that, with synergistic interactions between the spreading dynamics of two behaviors, both and can affect and the associated phase transition characteristically. To further demonstrate the role of the synergistic interactions, we show in Fig. 5 color coded values of in the parameter plane (, ) for , , , and . There are three regions in the parameter plane, determined by the two vertical lines at and , respectively, which are associated with characteristically distinct behavioral adoption dynamics. In region I (), only a small fraction of the nodes in layer adopt behavior 2. In region II (), there is a discontinuous phase transition, where a larger fraction of nodes adopt behavior for (white solid line). In region III (), there is a continuous phase transition. The distinct types of phase transition are predicted through our bifurcation analysis in Sec. III.

To gain further insights into the effects of synergistic interactions in behavioral adoption dynamics, we study the special case where the two types of behaviors are completely symmetric to each other. Fig. 6 (a) shows, for , , and , the dependence of on for different values of . In the absence of synergistic interactions, i.e., when the adoptions of behaviors and have no effect on each other, neither behavior can spread globally and either behavior can only be adopted by a small fraction of the nodes in the network. For (i.e., ), the nodes that have adopted behavior () only need additional pieces of information to adopt behavior (). As a result, the mutually cooperative spreading of behaviors and leads to a wide adoption of both behaviors. Increasing the synergistic interaction strength makes the dynamical correlation between the two layers stronger. The discontinuous phase is more clear when the network size is enlarged. However, the improvement in decreasing the deviation of the critical threshold is less, as shown in Fig. 6 (b). In this regime, the deviation is mainly because the theoretical method can not capture the strong dynamical correlation between layers.

### iv.2 General multiplex networks

We consider more general network topology for the network layers in the multiplex system, such as ER-ER, SF-SF, ER-SF and SF-ER. We use the standard configuration model Catanzaro et al. (2005) to construct SF networks with the degree distribution , where is the degree exponent and the coefficient is with the minimum degree and maximum degree . The average degrees of SF and ER networks are set as , and the network size is . For , e.g., and , we fix the final adoption size of behavior 1 and vary the type of network in layer .

To facilitate comparison, we set when layer is an ER network and if network is SF, so that the final adoption sizes of behavior 1 for both cases are approximately . As shown in Fig. 7(a), the network type in layer over which behavior 1 spreads has little effect on the spread of behavior 2. For the symmetric case , the dependence of on changes from being discontinuous to continuous as the network becomes more heterogeneous (i.e., SF) Wang et al. (2015a), as a strong heterogeneity makes it harder for nodes in the subcritical state to adopt a behavior simultaneously. Regardless of the network type, in general synergistic interactions can facilitate adoption of both behaviors and alter the nature of the associated phase transition.

## V Discussion

To understand social contagions in the human society at a quantitative level is of great importance in the modern time. While the spread of a single contagion can be analyzed through the traditional models of network spreading dynamics, the simultaneous presence and spreading of two or more contagions poses a challenge due to the mutual interplay between the underlying dynamical processes. As an initial effort to address this problem, we articulate a spreading model of multiple social behaviors on multiplex networks subject to synergistic interactions. For simplicity, we consider two-layer coupled networks and limit the number of distinct behaviors to two: one on each layer. The manifestation of the synergistic mechanism is that the adoption of the behavior by a node in one layer will increase the chance for the node that is simultaneously present in the other layer to adopt the behavior that spreads in that layer. The concrete setting enables us to develop an edge-based compartmental theory and a bifurcation analysis to uncover and explain how the synergistic interactions affects the spreading dynamics in terms of the final adoption size and the distinct phase transitions.

There are two types of synergistic interactions: asymmetric and symmetric. In the asymmetric case, the adoption threshold of one behavior in one network layer is less than that of the other behavior in the other layer. In this case, the adoption of the behavior with the higher threshold has no effect on the adoption of the other behavior. However, synergistic interactions can promote the adoption of both behaviors. In fact, the interaction strength and the information transmission rate of the behavior with the smaller threshold value can affect the nature of the phase transition of the behavior with the larger threshold: a small (large) value of the transmission rate of the former can lead to a discontinuous (continuous), first-(second-) order phase transition in the latter. In addition, a two stage spreading process arises: nodes adopting the small threshold behavior in one layer are more likely to adopt the large threshold behavior in the other layer, which stimulates the remaining nodes in this layer to quickly adopt the behavior. In the case of symmetric synergistic interactions, the adoption processes in both layers can affect each other on an equal footing. In this case, the interactions will greatly enhance the spreading of both behaviors in their respective layers through a first-order phase transition.

Many issues remain, such as the effect of heterogeneity in the synergistic strengths of the individual nodes on behavioral spreading and the impacts of degree correlation between the network layers. In general, there are two kinds of dynamical correlation: intralayer and interlayer. In each layer, the correlation can be described by the edge-based compartmental theory. To make a theoretical analysis feasible, we have neglected interlayer correlation, i.e., the dynamical correlation among nodes in distinct layers. However, in real situations, dynamical correlation may exist between the same node in different layers, depending on the strength of the synergistic interaction. If the interaction strength is not too large, interlayer dynamical correlation is weak. In this case, there is a good agreement between the theoretical prediction and the simulation results (e.g., Figs. 1 and 4). For relatively strong synergistic interaction (e.g., Fig. 6 for ), the simulation results deviate from the theoretical prediction. Increasing the size of network will not help reduce the deviation, as interlayer correlation can no longer be regarded as insignificant. A more accurate theory incorporating interlayer correlation is thus needed for synergistic affected information spreading in the strong interaction regime Wang et al. (2017).

## Acknowledgements

We are very grateful for the comments of anonymous reviewers. We thank A. Vespignani and Q. Zhang at the Laboratory for Modeling of Biological Socio-technical Systems (MOBS LAB) for valuable discussions and comments. This work was supported by the National Natural Science Foundation of China under Grants Nos. 11575041 and 61673086, the program of China Scholarships Council (No. 201606070059), and the Fundamental Research Funds for the Central Universities (Grant No. ZYGX2015J153). YCL would like to acknowledge support from the Vannevar Bush Faculty Fellowship program sponsored by the Basic Research Office of the Assistant Secretary of Defense for Research and Engineering and funded by the Office of Naval Research through Grant No. N00014-16-1-2828.

## References

- Barrat et al. (2008) A. Barrat, M. Barthélemy, and A. Vespignani, Dynamical Processes on Complex Networks (Cambridge University Press, Cambridge UK, 2008).
- Vespignani (2012) A. Vespignani, Nat. Phy. 8, 32 (2012).
- Pastor-Satorras et al. (2015) R. Pastor-Satorras, C. Castellano, P. Van Mieghem, and A. Vespignani, Rev. Mod. Phys. 87, 925 (2015).
- Valente (1996) T. W. Valente, Soc. Net. 18, 69 (1996).
- Christakis and Fowler (2007) N. A. Christakis and J. H. Fowler, New Eng. J. Med. 357, 370 (2007).
- Young (2009) H. P. Young, Ame. Econ. Rev. 99, 1899 (2009).
- Rogers (2010) E. M. Rogers, Diffusion of Innovations (Simon and Schuster, 2010).
- Centola (2011) D. Centola, Science 334, 1269 (2011).
- Banerjee et al. (2013) A. Banerjee, A. G. Chandrasekhar, E. Duflo, and M. O. Jackson, Science 341, 1236498 (2013).
- Zha et al. (2016) Y. Zha, T. Zhou, and C. Zhou, Proc. Nat. Acad. Sci. (USA) 113, 14627 (2016).
- Zhang et al. (2016) Z.-K. Zhang, C. Liu, X.-X. Zhan, X. Lu, C.-X. Zhang, and Y.-C. Zhang, Phys. Rep. 651, 1 (2016).
- Castellano et al. (2009) C. Castellano, S. Fortunato, and V. Loreto, Rev. Mod. Phys. 81, 591 (2009).
- Dodds and Watts (2004a) P. S. Dodds and D. J. Watts, Phys. Rev. Lett. 92, 218701 (2004a).
- Dodds and Watts (2005) P. S. Dodds and D. J. Watts, J. Theo. Biol. 232, 587 (2005).
- Granovetter (1973) M. S. Granovetter, Ame. J. Socio. 78, 1360 (1973).
- Watts (2002) D. J. Watts, Proc. Nat. Acad. Sci. (USA) 99, 5766 (2002).
- Gleeson and Cahalane (2007) J. P. Gleeson and D. J. Cahalane, Phys. Rev. E 75, 056103 (2007).
- Whitney (2010) D. E. Whitney, Phys. Rev. E 82, 066110 (2010).
- Hackett et al. (2011) A. Hackett, S. Melnik, and J. P. Gleeson, Phys. Rev. E 83, 056107 (2011).
- Zhuang et al. (2017) Y. Zhuang, A. Arenas, and O. Yağan, Phys. Rev. E 95, 012312 (2017).
- Gleeson (2008) J. P. Gleeson, Phys. Rev. E 77, 046117 (2008).
- Nematzadeh et al. (2014) A. Nematzadeh, E. Ferrara, A. Flammini, and Y.-Y. Ahn, Phys. Rev. Lett. 113, 088701 (2014).
- Brummitt et al. (2012) C. D. Brummitt, K.-M. Lee, and K.-I. Goh, Phys. Rev. E 85, 045102 (2012).
- Yağan and Gligor (2012) O. Yağan and V. Gligor, Phys. Rev. E 86, 036103 (2012).
- Cozzo et al. (2013) E. Cozzo, R. A. Banos, S. Meloni, and Y. Moreno, Phys. Rev. E 88, 050801 (2013).
- Dodds and Watts (2004b) P. S. Dodds and D. J. Watts, Phys. Rev. Lett. 92, 218701 (2004b).
- Krapivsky et al. (2011) P. L. Krapivsky, S. Redner, and D. Volovik, J. Stat. Mech. Theo. Exp. 2011, P12003 (2011).
- Zheng et al. (2013) M. Zheng, L. Lü, and M. Zhao, Phys. Rev. E 88, 012818 (2013).
- Liu et al. (2016) Q.-H. Liu, W. Wang, M. Tang, and H.-F. Zhang, Sci. Rep. 6, 25617 (2016).
- Wang et al. (2015a) W. Wang, M. Tang, H.-F. Zhang, and Y.-C. Lai, Phys. Rev. E 92, 012820 (2015a).
- Wang et al. (2016) W. Wang, M. Tang, P. Shu, and Z. Wang, New J. Phys. 18, 013029 (2016).
- Wang et al. (2015b) W. Wang, P. Shu, Y.-X. Zhu, M. Tang, and Y.-C. Zhang, Chaos 25, 103102 (2015b).
- Liu et al. (2017a) M.-X. Liu, W. Wang, Y. Liu, M. Tang, S.-M. Cai, and H.-F. Zhang, Phys. Rev. E 95, 052306 (2017a).
- Liu et al. (2017b) Q.-H. Liu, W. Wang, M. Tang, T. Zhou, and Y.-C. Lai, Phys. Rev. E 95, 042320 (2017b).
- Myers and Leskovec (2012) S. A. Myers and J. Leskovec, in Data Mining (ICDM), 2012 IEEE 12th International Conference on (IEEE, 2012) pp. 539–548.
- Zarezade et al. (2015) A. Zarezade, A. Khodadadi, M. Farajtabar, H. R. Rabiee, and H. Zha, arXiv preprint arXiv:1510.00936 (2015).
- Taubenberger and Morens (2006) J. K. Taubenberger and D. M. Morens, Rev Biomed 17, 69 (2006).
- Brundage and Shanks (2008) J. F. Brundage and G. Shanks, Emerg. Infec. Dise. 14, 1193 (2008).
- Chen et al. (2013) L. Chen, F. Ghanbarnejad, W. Cai, and P. Grassberger, EPL (Europhysics Letters) 104, 50001 (2013).
- Cai et al. (2015) W. Cai, L. Chen, F. Ghanbarnejad, and P. Grassberger, Nat. Phys. 11, 936 (2015).
- Hébert-Dufresne and Althouse (2015) L. Hébert-Dufresne and B. M. Althouse, Proc. Nat. Acad. Sci. (USA) 112, 10551 (2015).
- Domenico et al. (2016) M. D. Domenico, C. Granell, M. A. Porter, and A. Arenas, Nat. Phys. 12, 901 (2016).
- Kivelä et al. (2014) M. Kivelä, A. Arenas, M. Barthelemy, J. P. Gleeson, Y. Moreno, and M. A. Porter, J. Complex Net. 2, 203 (2014).
- Boccaletti et al. (2014) S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Genio, J. Gómez-Gardenes, M. Romance, I. Sendina-Nadal, Z. Wang, and M. Zanin, Phys. Rep. 544, 1 (2014).
- Catanzaro et al. (2005) M. Catanzaro, M. Boguñá, and R. Pastor-Satorras, Phys. Rev. E 71, 027103 (2005).
- Miller et al. (2012) J. C. Miller, A. C. Slim, and E. M. Volz, J. Roy. Soc. Interface 9, 890 (2012).
- Yang and Zhou (2012) Z. Yang and T. Zhou, Phys. Rev. E 85, 056106 (2012).
- Karrer and Newman (2010) B. Karrer and M. E. Newman, Phys. Rev. E 82, 016101 (2010).
- Baxter et al. (2010) G. J. Baxter, S. N. Dorogovtsev, A. V. Goltsev, and J. F. Mendes, Phys. Rev. E 82, 011103 (2010).
- Erdös and Rényi (1959) P. Erdös and A. Rényi, Pub. Math. (Debrecen) 6, 290 (1959).
- Barabási and Albert (1999) A.-L. Barabási and R. Albert, Science 286, 509 (1999).
- Wang et al. (2017) W. Wang, M. Tang, H. E. Stanley, and L. A. Braunstein, Rep. Prog. Phys. 80, 036603 (2017).