Modeling the Spread of Multiple Concurrent Contagions on Networks

Many contagions spread over various types of communication networks and their spreading dynamics have been extensively studied in the literature. Here we propose a general model for the concurrent spread of an arbitrary number of contagions in complex networks. The model is stochastic and runs in discrete time, and includes two widely used mechanisms by which a node can change its state. The first, termed the spontaneous state change mechanism, describes spontaneous transition to another state, while the second, termed the contact-induced state change mechanism, describes acquiring other contagions due to contact with the neighbors. We consider reactive discrete-time spreading processes of multiple concurrent contagions where time steps are of finite size without neglecting the possibility of multiple infecting events in a single time step. An essential element for making the model numerically tractable is the use of an approximation for the probability that a node transits to a specific state given any set of neighboring states. Different transmission probabilities may be present between each pair of states. We also derive corresponding continuous–time equations that are simple and intuitive. The model includes many well-known epidemic and rumor spreading models as a special case and it naturally captures spreading processes in multiplex networks.


Introduction
Epidemiological models, developed as tools for analyzing the spread and control of infectious diseases, have also been adapted across many scientific fields such as ecology, immunology, social science, computer science, marketing and economy. They focus on modeling the dynamics of contagious entities (also called ''memes'' in the literature) as diverse as communicable diseases, cultural characteristics (such as religious beliefs, fads or innovations), addictions, or information spread (through rumors, e-mail messages, web blogs, peer-to-peer computer networks, etc). Both deterministic and stochastic epidemic models have been suggested, addressing complementary questions [1][2][3][4]. Some well-known classical models are deterministic, and include, for example, the SIR (susceptible-infected-recovered) differential equation model of Kermack and McKendrick, which has proven useful in ascertaining gross factors affecting the rate of growth and the final size of an epidemic [5]. Stochastic models are preferable when studying a small community where the contact structure in the community contains small complete graphs, with households and other local social networks being common examples. But even when considering large communities, at which deterministic models primarily aim, some additional questions have been raised that can only be addressed with stochastic models [4]. Deterministic counterparts, working with the expected values of the corresponding stochastic models, are proposed for many of them to answer one of the fundamental questions for the propagation of a single contagion: will it infect a significant portion of the network or will it die out fast? Specifically, the existence of threshold values for the model parameters over which epidemic proportions can be reached has been studied. Earlier approaches have used a meanfield approximation [6], assuming homogeneous environments which are suited for simple network topologies. Afterwards, a heterogeneous mean-field (HMF) approximation has been introduced, which assumes that nodes with same degree behave in the same manner [7,8]. This approach can be applied to power-law networks, but its main assumption is not empirically or phenomenologically justified and it can result in different levels of accuracy [9]. One of the recent approaches, the so-called nonlinear dynamical system (NLDS) approach, describing the evolution of the probabilities of infection for every node [10][11][12][13][14] is widely accepted. As the number of states in the Markov chain which describes the dynamics of the whole network grows exponentially with the number of nodes, independence between the marginal distributions of the nodes is assumed in order to reduce the complexity of the models. This turns out to be a valid assumption in the vast majority of complex networks because the inherent topological disorder makes dynamical correlations not persistent.
Regarding the spread of more than one contagious entity, both deterministic as well as several stochastic models have been suggested [15][16][17][18][19][20][21][22][23][24][25][26][27][28]. However, the spreading rules in these models are specific to the problems the models are addressing. In [24] a SIR-like consecutive spreading of multiple viruses on special random networks has been introduced. It has afterwards been adapted for concurrent spreading of multiple viruses, and percolation analysis, which is suitable for SIR-like models, has been used to predict the epidemic sizes [25]. A generalization of the SIR model for two infections, the S Ã I 2 V Ã model, is proposed in [20], where a node's state is classified in one of the three categories: susceptible, infected or recovered (vigilant/vaccinated). The SI 1 I 2 S model is introduced in [20], where a susceptible node can become infected with or recover from one of two infections, and an extension is later made where it is possible that a node possesses both infections at the same time [19]. Spreading of multiple contagions has also been studied in the context of multiplex networks. Such examples are discussed in [27], [28] and [29], where the interplay between two different SIS propagations on two distinct layers of the multiplex network is considered.
Nonetheless, modeling the spread of multiple concurrent contagions with a discrete-time reactive process presents a significant distinct problem, and that is the question of how to deal with multiple simultaneous infecting events by different contagions in a single time step. Most of the proposed methods avoid this problem by assuming infinitesimally small time step sizes or asynchronous infections, allowing only a single event in a given time step. We are not aware of any work that studies the interplay between multiple competing contagions in networks and addresses the aforementioned problem directly, without neglecting the possibility of multiple simultaneous infecting events. The purpose of this paper is to propose such a model. Regarding the case with infinitesimally small time steps, we show that the differential equations of our model, derived from the discrete-time equations, naturally become simplified due to the lack of strong competition between the contagions on the level of a single node. Also, our work is among the first to propose a discrete-time stochastic model for the spread of an arbitrary, but finite, number of contagions which is general enough to describe a large class of spreading processes. The model we propose has two mechanisms of transition between two states of a node. The contact-induced transition mechanism is infection with other contagions due to contact with the neighbors, while the spontaneous transition mechanism characterizes a spontaneous transition to another model state without any contact with the neighbors. These two mechanisms encompass what is commonly met in the literature of modeling spreading processes and, as a result, our model generalizes some well-known single-and multiple-contagion spreading models.
As Daley and Kendall have pointed out, a mathematical model for the spreading of contagious entities can be (and has been) constructed in a number of different ways, depending on the mechanism postulated to describe the growth and decay of the actual spreading process [30]. Here we briefly describe two constructions commonly met in the literature. Consider a finite closed population of entities (agents) which is divided into three mutually exclusive and exhaustive classes x, y, z. If we assume that the only two transitions allowed are: from (x, y, z) to (x{1, yz1, z) at a rate proportional to xy, and from (x,y, z) to (x, y{1, zz1) at a rate proportional to y, then we obtain the deterministic Kermack-McKendrick epidemic model or stochastic SIR model, depending on what x, y, z are. On the other hand, assuming a transition from (x, y, z) to (x, y{1, zz1) at a rate proportional to yyzyz, one obtains the deterministic Daley-Kendall or Maki-Thompson model of rumor spreading. In the first construction, as is indicated by the transition rate xy, the growth in the number of entities is in one of the classes involved in the interaction. This particular interaction is between infected and susceptible individuals, and as a result, the number of infected individuals grows. In the second construction the increase in the number of entities in class z is additionally a result of interaction between classes unrelated to z, as is given with yy in the transition rate. In the Daley-Kendall and Maki-Thompson models this is based on the plausible hypothesis that an active spreader stops telling the rumor because when contacting another spreader it learns that the rumor has lost its news value. They both represent examples of the contact-induced transition mechanism. The spontaneous transition mechanism is also present in the first construction, where the transition rate from class y to z is proportional only to the number of members in class y.
The model suggested in this paper has three main characteristics. Firstly, it belongs to the class of stochastic discrete-time models, applies to arbitrary graphs, and can quantify the microscopic dynamics at the individual level by computing the probability that any given node is in a given state. A key instrument which we use to make this general model applicable for simulations is an approximation for the exact probability that a node will adopt a specific state from its neighbors. The approximation overcomes the challenge presented by the possibility of multiple simultaneous infections from the neighbors in a given time step which is a consequence of the finite sizes of the time steps.
Secondly, from the model one can derive its deterministic counterpart, both in difference and in differential equation form. Indeed, by assuming that the states of each node are independent random variables, one can derive a system of probability equations, which, in fact, represents a deterministic nonlinear dynamical system. Further, using a homogeneous or heterogeneous mean-field approximation for these deterministic dynamical systems, one can obtain models of differential equations describing the macroscopic spreading phenomena. Thirdly, the model generalizes the stochastic SIR, SIS and SIRS models for an arbitrary number of contagions or states and also suggests a stochastic microscopic Markov chain version of the deterministic Maki-Thompson model for an arbitrary number of rumors. It is general enough to capture spreading in multiplex networks as well.
We stress that although it is apparent that some, or even all of the real-world spreading phenomena do not fit into general simplified schemes and require special consideration of their details as they have characteristic modes of transmission, we believe that studying simple models may nevertheless be useful for understanding underlying principles of the spreading processes. At last, for the purpose of term unification, we shall refer to the contagious entities whose spread we model as infections or contagions, by analogy with the epidemiology literature, while also bearing in mind the generality of the model to describe the spreading of any other kind of entities through a network. We will often refer to them as states as well, from the representation of the dynamics of each node as a Markov chain.

Single-contagion spreading models
Before we define the model, we briefly turn our attention to the susceptible-infected-susceptible (SIS) model which we use as a paradigmatic model for the spread of one infection in a population of individuals connected in an arbitrary topology. The network of connections is represented by a simple, undirected and connected graph of N nodes whose adjacency matrix A~½a ij N|N is a binary valued matrix stating whether nodes i and j are connected (a ij~1 ) or not (a ij~0 ). An individual, which is represented by a node in the network, can be in either a susceptible (S) or infected (I) state. A susceptible node is healthy and it can receive the infection from infected neighbors. An infected node transmits the infection with probability b when contacting a neighbor. A successful transmission to a susceptible node causes it to become infected, and a successful transmission to an infected node has no effect. The probability that an infected node is cured and reverts back to the susceptible state is d.
Depending on the number of contacts a node makes, several spreading processes have been studied [12]. The one most often to be found in the literature is the reactive process. In this process an infected node contacts all of its neighbors and attempts to transmit the infection to each of them with probability b. Hence, a susceptible node can receive the infection from more than one neighbor. In this case, it chooses one of the successful transmissions and adopts the infection transmitted by that contact. Since the SIS model describes the spread of only one infection in the network, it is irrelevant which particular successful contact a node chooses in order to adopt the infection.
A discrete-time stochastic mathematical model of the thus described process is as follows. The state of node i at time t is described by a state vector containing a 1 in the component corresponding to the current state of the node, and 0 in the other: The probability mass function which states the probability of being in each of the states is given with the probability vector The model equations describing the time evolution of the probability vector for each node are: where f I i (t) is the probability that node i receives the infection from at least one infected neighbor. When writing an expression for f I i (t), it is commonly assumed that transmission events in the current or past time steps are independent of each other. The specific spreading process determines the form of f I i (t) and for the reactive process it reads The product in (1) is the probability of the event that none of the infected neighbors (a ij s j,I (t)~1) transmits the infection to node i. Hence, f I i (t) is in fact given with the probability of the opposite event, which is that at least one infected neighbor manages to transmit the infection to node i. The state diagram that summarizes the Markov chain diagram for the dynamics of a single node in the SIS model is given in Fig. 1.
In order to infer the analogy of this (reactive) spreading process where a single contagious entity spreads in a network to the one where more contagious entities spread simultaneously, we use the binomial theorem to rewrite (1) for a node i with d infected neighbors: Equation (2) clearly shows that the probability f I i (t) of receiving the infection from the neighbors is the sum of the probability that exactly k neighbors have successfully transmitted the infection to node i, where k goes from 1 of the infected neighbors to all d of the infected neighbors. The binomial coefficient takes into account all combinations of k successful transmissions out of possible d transmissions. The term (1{b) d which corresponds to the event that no infected neighbor successfully transmitted the infection is canceled out in (2).
Here we stress that when exactly k neighbors have successfully transmitted the infection, node i chooses only one of these, which means that the probability of receiving the infection from one of those k successful transmissions is 1 there are k successful transmissions, the probability of receiving the infection from any one of those is which is how the terms in the sum in (2) are obtained. This happens since all transmission events have the same probability b of occurring, which, in turn, leads to the possibility to write f I i (t) in the product form as in (1).
On the other hand, when more contagions spread in the network simultaneously and have different transmission probabilities, the fact that a node chooses only one of the k successfully transmitted, generally different, contagions makes it impossible to write the corresponding probability of infection f I i (t) in the product form as in (1). The expression for f I i (t), as we shall see in the following, will be analogous to (2), where all combinations of neighbors infected with different contagions will be included. On the other hand, infected nodes spontaneously recover with probability c, and become susceptible again. The contact-induced transition mechanism is represented by a curvy arrow, whereas the spontaneous transition mechanism is represented by a less curvy arrow. doi:10.1371/journal.pone.0095669.g001

Model description
We study the following discrete-time stochastic process for the spreading of multiple contagions in a network. The network is represented by a simple, undirected and connected graph G~(V ,E) with node set V and link set E. The graph's adjacency matrix is A~½a ij N|N , where N is the number of nodes in the network. We assume that every node i f1,2, . . . ,Ng in the network is in a certain state k f1,2, . . . ,mg. A state k can represent the healthy state or the recovered state, as in the classical epidemic models, or it can denote that the node is infected with a certain contagion k. Node i's state vector at time t is represented by where s i,k (t), k~1, . . . ,m, is a Bernoulli random variable indicating whether node i is in state k. A node can be in only one of the model states at a time, meaning that only one component of the state vector is 1 and all others are 0 for each t. This is known as the 1-of-m coding scheme. The probability mass function corresponding to the state vector s i (t) is where p i,k (t) is the probability that node i is in state k at time step t. Naturally, P m k~1 p i,k (t)~1 holds. The current state of a node can be changed by one of two mechanisms, as is illustrated in Fig. 2. The first mechanism is spontaneous transition to another state, without any contact with the neighbors. A node spontaneously abandons state k and transits to state l with probability d kl ½0,1. Note that P m l~1 d kl ƒ1 and that in general d kk =0. The process is equivalent to the rolling of a loaded (mz1)-sided dice, where the outcome of the roll is the state to which node i spontaneously transitions, and the (mz1)th side is the event that node i does not spontaneously change its state. This mechanism is analogous to node curing in epidemic models where nodes change from an infected to a healthy state. Should the node not make a spontaneous change, it proceeds with the second mechanism. The second mechanism is transition to other state due to contact with the neighbors and we also refer to it as the contact-induced state change mechanism. This mechanism is analogous to node infecting with a contagion from its neighbors. The transmission probability b kl [½0,1, is the probability that a node in state k will change its state when contacting a neighbor in state l. After a successful transmission, the receiving node may transit to state l or, stimulated by the communication, adopt another state. To simplify the modeling process, we make the usual assumption that transmission events between separate pairs of nodes are independent of each other.
The model equations are where p(e h ,c q )~P Equation (3) describes the time evolution of the probability that node i will be in state l in the next time step. The first term of the right-hand side of (3) gives the probability with which node i spontaneously transits to state l from its current state (first panel of Fig. 2). The second term encompasses the probability that, provided no spontaneous transition occurs, node i will adopt state l due to the second mechanism of state change (second and third panel of Fig. 2). The probability of this event is 1{ is the probability that node i which is currently in state k transits to state l due to contact with the neighbors. The third term states the probability with which a node currently in state l is not affected by any of the two mechanisms of state change. This event happens when the node does not spontaneously transit to any other state, and does not receive any other state or contagion when contacting the neighbors, the probability of which is g 0 i,l (t) (fourth panel of Fig. 2). The probability g l i,k (t) of state transmission, or infection with contagion l, due to communication with the neighbors, given with (4), takes into account the two possible constructions of the contact-induced state change mechanism described previously in the introduction. The first one is typical for well-known epidemic models where node i can transit to state l directly by contacting a neighbor in state l, and this is depicted in the second panel of Fig. 2. The second one is encountered in social spreading processes, where the exposure of node i to another state h upon contact with the neighbors stimulates it to adopt state l. This is illustrated in the third panel of Fig. 2. Such is the case, for example, in the Maki-Thompson model where, upon contact of two informed nodes, one becomes a stifler since it loses interest in the rumor it possesses. Another, more complex, example would be the case where a single state is represented as a set of multiple distinct features, like in the language evolution models in [31,32]. Upon successful infection from neighbor j, node i may copy only a fraction of the features from node j that it does not possess itself, effectively adopting a state that differs from node j's state. A similar example is presented in [19], where exposure of an infected node to the other infection in the network causes it to transit to a third state which signifies that it possesses both infections simultaneously.
is the probability that state h will cause change to node i in state k from any combination of its infectious neighboring links and is the analogue of (1) for the SIS model. t l kh [f0,1g is an indicator variable which states whether the exposure of a node in state k to state h will stimulate it to adopt state l. Note that when node i in state k transits to state l only by receiving it directly from its neighbors, t l kh~1 for h~l and t l kh~0 for all h=l. Hence, we have g l i,k (t)~f l i,k (t), which is the most common case in the existing models. The sum in (4) goes over all states h which can stimulate a node in state k to adopt state l. In essence, g l i,k (t) is a generalization of f l i,k (t) for the purpose of encompassing the second construction. A slightly more general case is the one when t l kh are parameters such that P m l~1 t l kh~1 and t l kh [½0,1. This makes t l kh , for l~1, . . . ,m, to act as weights to the transitioning states instead of indicator variables. Again, a good example would be the case with states which describe the possession of multiple features, where upon successful infection a node may have multiple possible courses of action, regarding the choice of copying individual features. We do not treat such cases in this paper.
The specific spreading process determines the form of f l i,k (t). In this paper we concentrate on generalizing the reactive process as indicated in [12] for an arbitrary number of node states. As mentioned previously, this is the type of process most often to be found in the literature. In the reactive process, a node contacts all of its neighbors in every time step and tries to spread its current contagion to all of them, i.e. it tries to convince the neighbors to adopt the state that the node is currently in. As stated above, the probability of transmission b kl depends on the states k and l of the contacting nodes (receiving and sending node, respectively). The link over which a successful transmission has occurred is said to be infectious at the given time step. As transmissions are independent of each other, multiple infectious links may occur in a single time step at a receiving node. Most of the papers in the literature avoid this problem by assuming infinitesimally small time steps, virtually avoiding multiple simultaneous infection events. For node i with a Solid colored arrows indicate successful state transmissions, i.e. infectious links, and dashed lines indicate an unsuccessful state transmission. The probabilities of the realized transmission events are depicted next to each line. Solid gray lines indicate that the nodes have not been in contact at the given time step; a spontaneous transition has taken place instead. From top to bottom panel, descriptions go as follows. Panel 1: node i changes its state spontaneously to state 2 after previously having been in state 1. The probability of state change with this mechanism is d 12 . Panel 2: Node i does not make a spontaneous transition, and changes its state as a result of getting infected with state 3 from its neighbors. Note that a neighbor in state 2 also makes successful transmission, however, node i chooses state 3 transmitted from one of the other two successful transmissions. The probability of state change with this mechanism is (1{ . Panel 3: node i changes its state as a result of getting infected with state 4 from its neighbors, a contact which stimulates it to adopt state 2. The probability of state change with this mechanism is (1{ 2 14 , where f 4 i,1 (t)~(1{b 12 )(1{b 13 ) 2 b 14 and t 2 14~1 . Panel 4: node i maintains its state since none of the two mechanisms of state change caused it to make a transition. The probability of this event is the product of the probability that no spontaneous transition occurs and no state is transmitted upon contact with the neighbors (1{  Now, the probability f l i,k (t) that node i which is in state k adopts state l from any combination of its infected neighbors for the reactive process is given with (5). Equation (5) goes over every possible configuration of states c q at the neighbors of node i and every possible event e h . The terms s cq j,l (t) and ½e h j are used to filter only the combinations where state l is involved. s cq j,l (t) is simply the Bernoulli random variable which indicates whether node j is in state l in the configuration c q . When there are multiple successful infectious links in e h , i.e. a total of P di j~1 ½e h j infectious links, node i chooses one of the transmitted states to adopt, with uniform probability In order to provide the expression (6) for p(e h ,c q ) we denote the current state of each neighbor j with r j (t). Since we make the assumption that a transmission between contacting nodes is independent of other transmissions in the network in the current or past time steps, the probability that event e h occurs at node i is simply a product of the probabilities that the respective links in the event have become infectious or not. b krj is the probability that the link with neighbor j becomes infectious. The argument from r j (t) is omitted for brevity (we make the same omissions in the rest of the text).
Lastly, g 0 i,l (t) in (3), given with (7), is the probability that node i in state l does not adopt any state from its neighbors, i.e. no infectious links occurred in the current time step.
The computational complexity of (5) is of order d i (2m) di . Due to the exponential dependence on the degree d i of node i, the probability f l i,k (t) becomes numerically intractable for nodes of high degree. Therefore, further in the text we give an approximation for (5) which allows the application of the model for networks with high-degree nodes. We also speculate that the difficulty in finding a suitable expression for the probability f l i,k (t) which is numerically feasible is the reason why a general model of this kind has not appeared so far.
The model has three types of parameters: the transmission probabilities B~½b kl m|m , the probabilities of spontaneous transition D~½d kl m|m , and the indicator variables T~½t l kh m|m|m . Although the number of parameters is large, for most of the models in the literature that we generalize, these matrices adopt a sparse structure. In the rest of the paper we continue with the implicit assumption we have made so far, that the parameters are the same for every node. However, node dependence of the parameter values is easily incorporated by making the parameters different for each node.

Approximation
The complete derivation of the approximation, which is given in Section S1 of Text S1, leads to the following approximate expression for the probability f l i,k (t) given with Eq. (5): The derivation process was conducted while bearing in mind the compatibility with the deterministic counterpart, where instead of state vectors there are probability vectors. Hence, Eq. (8) can naturally be used for the deterministic case. As a result, it takes slightly complex form. However, we can make a simplification for the stochastic case. Let d Observe that when m kh~mkq , i.e. b kh~bkq , for all h,q, we havẽ f f l i,k~f l i,k for all l. We stress that the term which represents the fraction of neighbors in state l, is primarily the term that is approximated. As already mentioned, for models where only one state is being transmitted by the contact-induced mechanism, this fraction is equal to 1 and the probability of receiving the state can also be written in the product form (1). The nonlinearities in the model arise from the product term in (6). In existing models for the spread of a single contagion, this product is most often linearized using a general form of the Weierstrass product inequality for the purpose of model analysis. Particularly, Eq. (1) is usually substituted by b X N j~1 a ij s j,I (t), an approximation which holds only for b%1.
An assessment of the accuracy of the approximation is presented in the section where we discuss our generalization of the SI 1 I 2 S model and further in the section where we present an example numerical model. We compare the non-approximated and approximated version of the deterministic counterpart of our model, described in the following section. Additional assessments are also given in the supporting material (Section S2 of Text S1). Results indicate that the absolute value of the error is of order 10 {4 to 10 {2 for a single node, depending on the specific scenario. Increasing the number of neighbors with one of the states produces a lower precision of the approximation when the state transmission probabilities are high (w0:25) than when they are low. Also, the approximation usually overestimates the actual probabilities for states with high transmission probabilities when they interplay with states with low transmission probabilities, for which the actual probabilities are underestimated as a result. On the other hand, the approximation error is very low when the set of state transmission probability values has a small variance. An important result is that the error does not seem to accumulate over time and individual errors do not influence each other significantly. This has been observed in the SI 1 I 2 S model and its altered version, as the predictions of the fixed points produced by the approximated version are in agreement with those produced by the non-approximated version. Further, the numerical example showed that the approximation is accurate throughout the time evolution even when the spreading dynamics showed oscillatory behavior. This was also the case for a network with high degree nodes.
Besides making the model numerically tractable, the approximation also allows for the estimation of the model parameters and their joint posterior distribution when data about the process being modeled are available. The joint posterior distribution of the parameters is a high-dimensional distribution which can be approximated by Markov Chain Monte Carlo sampling methods, such as the the Metropolis-Hastings algorithm. The posterior distribution is proportional to the likelihood of the observed data (in this case, the states of the nodes in the network at certain time steps), and the probability mass vectors for each node (obtained by running Eq. (10)) are required to calculate the likelihood. The probability mass vectors, in turn, depend on the probabilities f l i,k (t) which can easily be calculated using the approximationf f l i,k (t). Using Bayesian inference, an analysis of any function of the parameters can be performed.

The deterministic counterpart of the model and its continuous-time form
By applying the expectation operator to Eq. (3), and taking into account the assumption of independence of transmission events in the current or past time steps, as well as the fact that p i,k (t)~E½s i,k (t), from Eqs. (3)-(7) a deterministic discrete-time version of the model is obtained, which is basically a discrete-time nonlinear dynamical system: where For brevity, we only presented the approximated equations of the model here, while the non-approximated version is analogous to Eq. (5). This deterministic discrete-time version of the model describes the dynamics of expected-value quantities of the stochastic model. Instead of running just one stochastic realization with the stochastic version of the model, or running sufficiently many to produce average results, one can use the deterministic model to obtain the average dynamics of the spreading process and make predictions as to its future state. Hence, we expect that the deterministic results will be comparable with those obtained from the mean values of the stochastic realizations in the limit of large network sizes. Furthermore, it also allows for a simpler analysis of the dynamical behavior of the model, which we leave for future work. For example, using a classical result for the weak ergodicity of time-inhomogeneous Markov chains by Wolfowitz [33], one can determine conditions by which the deterministic version of the model has a globally stable fixed point, which means that the average dynamics of the stochastic model stabilize. Such an analysis is given in [20], for instance.
In order to obtain a deterministic continuous-time, i.e. differential equation, model, the product in (13) is typically linearized, which holds for b lh %1,V l,h. Rearranging the terms in Eq. (10) in a way that one can calculate the limit dp i,l (t) dt~l im Dt?0 p i,l (tzDt){p i,l (t) Dt ,Vl, the deterministic differential equations for the evolution of the probability mass function of node i are obtained (derived in detail in Section S3 of Text S1): or in the case when t h kh~1 , Vh, it reads where each parameter that represents a probability in the discretetime model is substituted by its respective rate (e.g.b b kl Dt~b kl , for time step Dt). The first and the second term on the right hand side of Eqs. (14) and (15), which increase the rate of change of p i,l , correspond to transitions from other states to state l due to the spontaneous and contact-induced state change mechanism respectively, while the other two terms correspond to the respective state transitions in the opposite direction. One can also note that the second term on the right hand side of Eq. (15) has a simpler form than the corresponding term from the discrete-time case (Eqs. (10)- (13)). Specifically, the continuous-time one depends only on the components of the neighbors' probability vectors that correspond to state l. This is due to the infinitesimally small time step size which prevents the occurrence of more than one transmission event in any given time instance. Equation (15) written in matrix form reads where the operator 0 represents the Hadamard (or element-wise) product. The matrix P has the probability vectors of the nodes as rows, while the matricesD andB B have the rates of the spontaneous and the contact-induced transition mechanism as elements, respectively. 1 N|m is an all-ones matrix with N rows and m columns. A is the adjacency matrix of the network, as mentioned earlier.

Special Cases of the Model
In this section we present some well-known models which the proposed model generalizes. Since the model equations (3) constitute an inhomogeneous Markov chain for each node, we give the state diagrams of the Markov chains for each of the models in Figs. 1-6. The models discussed are the widely known SIS, SIR and SIRS epidemic spreading models; the Maki-Thompson rumor spreading model; and the SI 1 I 2 S model where two contagions spread concurrently in the network, for which we also assess the accuracy of our approximation. We present the natural extension of our model on multiplex networks, as well. Although the model equations appear to be complex in the general case, we can observe that they reduce to much simpler ones due to the sparsity of the parameter matrices of the presented models. For the discrete-time forms of each of the models, the parameters which represent probabilities are used in Eq. (3) and Eq. (10), while for the continuous-time forms the respective rate parameters are used in Eq. (14). Other models which we have met in the literature and which the proposed model generalizes are given in [19,20].

The epidemic spreading models
The SIS model. The SIS model is described in the previous section. Its state diagram is presented in Fig. 1. Recall that the state of node i at time t is described by a state vector of length m~2 for the SIS model, specifying whether it is in the susceptible (S) or infected state (I). The model equations (3) for the evolution of the probability mass function become:  On the other hand, infected nodes spontaneously recover with probability c. However, they only obtain temporal immunity which they lose with probability a. The contact-induced transition mechanism is represented by a curvy arrow, whereas the spontaneous transition mechanism is represented by a less curvy arrow. doi:10.1371/journal.pone.0095669.g004 An ignorant node can become a spreader by contacting its neighbors that spread the rumor, with rate b. On the other hand, spreader nodes become stiflers by contacting other spreaders or stiflers with rate a, or by spontaneously transitioning to the stifler state with rate d. The contact-induced transition mechanism is represented by a curvy arrow, whereas the spontaneous transition mechanism is represented by a less curvy arrow. The dashed arrow denotes that t R SS , i.e. spreader becomes stifler by contacting other spreaders with rate a. doi:10.1371/journal.pone.0095669.g005 Figure 6. State diagram of a node for the SI 1 I 2 S model. The diagram shows the dynamics of a single node. Curvy arrows depict state change due to contact with the neighbors, while less curvy arrows depict spontaneous state change. doi:10.1371/journal.pone.0095669.g006 Now, since the contact-induced transition mechanism of the SIS model is infection spread from infected nodes to susceptible nodes, and the spontaneous transition mechanism involves spontaneous transition only from the I to the S state, the parameters of the model are: The parameters for the continuous-time case,B B andD, are defined analogously. Note that in the SIS model no interaction with a given state stimulates a node to adopt a state other than the given one involved in the interaction. Therefore, t h kh~1 , Vk,h, so we have t I SI~1 . This means that g I i,S (t)~f I i,S (t). Also, as noted previously, when a node can be affected by only one infection through the contact-induced mechanism, no approximation is needed to write out (5) and f I i,S (t)~f f I i,S (t) can be written in the product form as in (2). This gives g S i,I (t)~0, g S i,S (t)~0, g I i,I (t)~0, and g I i,S (t)~f f I i,S (t)~1{g 0 i,S (t)~1{ P N j~1 (1{ba ij s j,I (t)). Taking all of the aforementioned into account, we find that (17) reduces to the known equations for the SIS model: Similarly, using the parametersB B andD with Eq. (14), we obtain the continuous-time equations of the SIS model: The SIR model. The SIR, or susceptible-infected-recovered model, is used to describe the spread of an infection for which permanent immunity is obtained after the end of the infectious period. The state diagram that describes its Markov chain is given in Fig. 3. A node can be in one of three states: the susceptible (S), infected (I) and recovered (R) state, hence s i (t)~s i,S (t) s i,I (t) s i,R (t) ½ for a given node i. The general model equations (3) for the SIR model are: According to the description of the model, the parameters of the model are Similarly, using the parametersB B andD D with Eq. (14), we obtain the continuous-time equations of the SIR model: The SIRS model. The SIRS, or susceptible-infected-recovered-susceptible model, is used to describe the spread of an infection for which temporary immunity is obtained after the end of the infectious period. The state diagram that describes its Markov chain is given in Fig. 4. It is very similar to the SIR model, the only difference being the addition of a spontaneous transition link from the recovered to the susceptible state that describes the temporary immunity. Having the same states as the SIR model, the general model equations (3) for the SIRS model coincide with those of the SIR model, Eq. (20). As mentioned earlier, the only difference in the model parameters is that d RS~a in the SIRS model, whereas d RS~0 in the SIR model. Hence, the equations (3) reduce to p i,S (tzDt)~s i,R (t)azs i,S (t)P N j~1 (1{ba ij s j,I (t)) p i,I (tzDt)~s i,I (t)(1{c)zs i,S (t)½1{ j~1 (1{ba ij s j,I (t)) p i,R (tzDt)~s i, Similarly, the continuous time equations of the SIRS model read The Maki-Thompson model for rumor spreading The Maki-Thompson model [34] is a popular variant of the classical model of rumor spreading of Daley and Kendall [35]. In both models, which operate in continuous time, each node can be in one of three different states: ignorant (I), spreader (S) and stifler (R). Nodes in state I are uninformed, or ignorant, of the given rumor and are thus susceptible to it. Nodes in state S actively spread the rumor, while nodes in state R are aware of the rumor, but they have lost interest in it and no longer spread it. When an ignorant node contacts a spreader, it becomes a spreader as well with rate l. If, on the other hand, a spreader contacts a stifler or another spreader, it becomes a stifler at a rate a. The two models, Maki-Thompson and Daley-Kendall, differ in the contact mechanism. The Daley-Kendall model adopts a pair-wise contact mechanism, i.e. for a given neighboring pair only one communication is assumed, which can affect both nodes simultaneously. For example, two neighboring spreaders either both become stiflers due to a single successful transmission or otherwise they both remain spreaders. Our model does not generalize such models. On the other hand, the Maki-Thompson model adopts a directed contact mechanism, i.e. only the node which initiates the contact can change its state. Hence, separate contacts will be initiated for both directions in every time step for the discrete time case, as we work with a reactive process. Both models assume a homogeneously mixed population or an undirected network; however, the Maki-Thompson model can easily be applied on directed networks as well. Also note the resemblance between the ignorant, spreader and stifler state in the Maki-Thompson model and the susceptible, infected and recovered state in the SIR model. The difference with the epidemic models is the part of the contactinduced transition mechanism which comes from contacts initiated by the spreaders. The mechanism of spontaneous transition to another state is absent from the original versions of both models.
A more recent study has supplemented both models with a mechanism of spontaneous rumor forgetting with rate d by which spreader nodes can also become stiflers [36]. The state diagram describing the Markov chain of this version of the Maki-Thompson model is given in Fig. 5. Here we show that the proposed model reduces to the Maki-Thompson model with rumor forgetting. Moreover, we recover the differential equations that have actually been used in [36]. The state of a node is described by a vector of length m~3: s i (t)~s i,I (t) s i,S (t) s i,R (t) ½ , and so is the probability mass function which gives the probability of being in each state. Regarding the model parameters, note that l, a and d are rates, and the corresponding probabilities are obtained by multiplying these with the length of the time step Dt in which a contact occurs. Hence, the parameters are  This means that f I i, in general. The contact-induced transition mechanism in the Maki-Thompson (and Daley-Kendall) model describes the increase in the number of both the S-nodes and the R-nodes. The increase in the number of spreaders is achieved by the interaction of nodes in states I and S, which is an infectious spread of the S state. The increase in the number of stiflers occurs because of interactions between nodes in states S and R, which can be viewed as an infectious spread of the R state, and because of interactions between nodes in state S, where the contact stimulates a node to adopt the R state. Therefore, t h kh~1 , Vk,h, except for k~h~S, where t S SS~0 and t R SS~1 , which results in g S i,R (t)~0, (1{aDt a ij s j,S (t)): The probabilities g 0 i,I (t) and g 0 i,S (t) of not changing the current state given with (7) become g 0 i,I (t)~P N j~1 (1{lDt a ij s j,S (t)), 1{aDt a ij (s j,S (t)zs j,R (t)) Â Ã , and g 0 i,R (t)~1, which gives rise to the term s i,R (t) in the last equation of (25).
Regarding the continuous-time case, using the rate parameter matrices,B B andD D, with the matrix T and Eq. (14), we obtain the following deterministic differential equations for the evolution of the probability mass function of node i, which are the same as the ones from [36]: Further, if we use a homogeneous mean-field approximation, which means that we assume that every node i has the same degree k k~X N j~1 a ij , and hence the same dynamical behavior, i.e. p i (t)~p(t) for all i[f1, . . . ,Ng, one obtains dp I dt~{ l k kp I p S dp S dt~l Lastly, the original version of the Maki-Thompson model is obtained for d~0.
The SI 1 I 2 S model As a suitable example for assessing the accuracy of our approximation, we present our discrete-time generalization of the continuous-time SI 1 I 2 S model [26]. The state diagram that describes a node's transition probabilities is presented in Fig. 6. In this model, there is a competition between two different contagions (I 1 and I 2 , or states 1 and 2, respectively), in the sense of which contagion will infect a given susceptible node (in state S, or state 0). Their respective transmission probabilities are b 1 and b 2 . Also, once infected, a node can become susceptible again with a probability d 1 or d 2 , with respect to I 1 and I 2 . Taking all of this into account, the discrete-time model equations (10) for the deterministic SI 1 I 2 S are wherẽ for l[f1,2g, and We can derive the corresponding continuous-time model equations (15) as well, which read Equation (29) coincides with the one derived in [26].
In order to assess the accuracy of the approximation of our model for the discrete-time case (Eqs. (10)-(13)), we compare it with the actual non-approximated version on several small networks (complete, star, ring and lattice graphs). The limitation on the network size comes as a result of the computational complexity of the non-approximated version. We calculate the macroscopic fixed points, i.e. percentage of nodes in each of the states, produced by both versions of the model for a given set of parameters. The norm of the error vector for every combination of b 1 and b 2 is calculated. As the approximation refers to the contactinduced state change mechanism, we keep d 1 and d 2 fixed at 0:075 and 0:05, respectively. The results show that our approximation produces the same macroscopic fixed points as the non-approximated version, for almost every combination of b 1 and b 2 . This can be expected, since in this model, as shown in [26], there is always a clear winner except for the case when b 1 d 1~b 2 d 2 , which is the line where our approximation does not produce zero error. This is shown in Fig. 7, for the complete graph, the star graph and the lattice graph.
A model with more complex behavior can be easily constructed, and it will help in assessing the accuracy of our approximation. We alter the aforementioned SI 1 I 2 S model by replacing the spontaneous transition links with contact-based transition links. The state diagram of this model is presented in Fig. 8. The process that we model can be thought of as a competition between two political parties (I 1 and I 2 ). As before, a person that does not support either party can become a supporter by contacting other supporters of one of the parties. Also, a voter can decide not to vote by contacting other such persons, having previously supported one of the parties (unlike the previous model, where such state change can only be made spontaneously). The macroscopic fixed points were compared for the approximated and the non-approximated version of our model, for every combination of b 01 and b 02 . Similarly as before, we kept b 10 and b 20 fixed at 0:075 and 0:05, respectively. In this model a clear winner is not always found. The approximation refers to the contact-induced state change mechanism of the susceptible node, the same as in the original SI 1 I 2 S model. The results of the comparison between both versions, presented in Fig. 9, demonstrate how much the individual approximation errors influence each other and accumulate over time. They clearly indicate the regions where our approximation works best. For relatively comparable values as well as for small values (v0:25) of both of the parameters b 01 and b 02 our approximation produces a negligibly small error on every graph. We find these results satisfactory as working within these regions has been standard practice in the literature. The error in the other regions is low as well. The results for the ring graph are not shown, since the error is significantly smaller than the errors for the other graphs.

An example model for the spread of three innovations
While the previous example model was suitable for presenting the regions where our approximation very closely matched the non-approximated deterministic version of our model, the macroscopic behaviors always converged to a fixed point where one of the states had vanished from the system. We observed that the deterministic behaviors did not show oscillations, due to the simplicity of the model. In order to test the approximation on a system with oscillatory behavior, we introduce an example model for the spread of three innovations, whose state diagram is shown in Fig. 10. Here, each node possesses one of the three innovations (there is no susceptible state in this model), and contact-induced  transitions can occur between any two states, including transitions to the same state. The parameter matrix B that was used is Relatively large values were selected to establish whether our approximation works well even for the parameters where the error was not zero in the previous assessments (Fig. 9). Also, as mentioned earlier, they were selected to create oscillatory, i.e. cyclic-like, behavior. More specifically, should we leave only the maximal element of each row of the matrix B, the respective state diagram will depict a cycle graph. The link widths in the state diagram ( Fig. 10) are proportional to the corresponding parameters from the matrix B as well. A sample execution on a lattice network with periodic boundary conditions of 65536 (256|256) nodes is presented in Fig. 11, that demonstrates the oscillatory behavior. Snapshots of the system state at eight different time steps are given.
Due to the large number of parameters in this model, a comparison of the fixed point values is inappropriate. However, another suitable way for assessing the accuracy of our approximation is to compare the actual simulations of the stochastic model, which are described by the non-approximated version of the model (Eqs. (3)-(7)), with the approximated version of the model (Eq. (9)). While comparing the non-approximated and approximated deterministic models is not possible for large networks and networks with high-degree nodes (due to the computational complexity of the non-approximated version of the model), it is easily done with the stochastic model. For a lattice network with periodic boundary conditions of 16384 (128|128) nodes, the comparison results are shown in Fig. 12, left panel. They were produced by averaging over 1000 executions. The markers display the execution results of the approximated version, while the lines display the results of the non-approximated version, i.e. the actual model simulations. Although we selected relatively large values for the parameters, our approximated version matched the non-approximated version very closely. To prove that this is also the case for networks that have nodes with large degrees, we performed the same comparison on a power grid network [37] of 4941 nodes and 13188 links, whose largest degree is 19. Results, shown in Fig. 12, right panel, again indicate that the approximation works well. At last, we performed similar comparisons on models with smaller parameter values, and the results, which we do not present here, showed that the approximated version matched the non-approximated version even better.

Spreading in multiplex networks
To demonstrate the broad applicability of our model, we present the multiplex adaption of our model as a special case. A multiplex, or composite network, is a network whose nodes are interconnected with various types of links (Fig. 13, first panel) [38]. Each link type corresponds to a different layer in the network (Fig. 13, second panel), however, a node can only be in one single state over all layers in a given time step. Depending on how the multiplex network is defined, different variations may exist, e.g. a node's state can be depicted by a separate state for each layer. We focus on the single state case. As before, multiple contagions or states can spread in the network. However, the transmission  probabilities between the states are layer-dependent, i.e. each of the contagions spreads differently on a different layer. Specifically, we have multiple adjacency matrices A f1g ,A f2g Á Á Á A fLg that describe the connectivity in each of the L layers of the multiplex network. We also have multiple contact-induced transition mechanism matrices B f1g ,B f2g Á Á Á B fLg , one for each layer, respectively, as shown in Fig. 13, third panel. Hence, the multiplex spreading equation for the probability that node i will be in state l aligns with Eq. (10) withf and Figure 11. Snapshots of one sample execution of the I 1 I 2 I 3 model on a lattice network with periodic boundary conditions of 65536 (256|256) nodes. The snapshots were taken at different time steps, as indicated below each individual snapshot. Cyclic-like behavior is clearly seen. doi:10.1371/journal.pone.0095669.g011 As can be seen, the model is powerful enough to capture spreading in multiplex networks. One example of such spreading is the multiplex variant of the SI 1 I 2 S model [27]. There, two contagions (or memes) spread on two separate layers of the multiplex network. Mutual exclusivity is assumed between the contagions, i.e. a node infected with one of the contagions is immune to the other. The state diagram describing the model is presented in Fig. 13, third panel, and it coincides with the previous example. However, in this case, an infection has zero transmission probability on the layers it does not spread on. Hence, b The corresponding differential equations read

Conclusions
In summary, this work has two main contributions. Firstly, we have proposed a general model for the spread of an arbitrary number of infections or contagions on networks, in which several contagions can simultaneously compete to infect a node. The model is stochastic and runs in discrete time, since difference equations are used to describe the evolution of the probability of being in each state. It can describe the spread of not only biological infections, but also social contagions such as information and rumors, cultural characteristics such as innovations and languages, and other systemic entities whose spread on networks can be described by the two mechanisms of state change incorporated in the model, which we refer to the spontaneous and the contact-induced state change mechanism. The first mechanism describes spontaneous transition to another state, and corresponds to the curing mechanism in classical epidemic models, while the second mechanism describes infection with other contagions or states due to contact with the neighbors, and corresponds to the spreading process in classical epidemic models. Secondly, an essential step for making the model applicable for simulations on networks is that we use the approximation (8) for the exact probability (5) that a node will adopt a specific state from its neighbors, which may possess any of the states in the network. The approximation showed high accuracy in most of the tested cases.
The proposed model generalizes classical epidemic models such as the SIS, SIR and SIRS models. Additionally, the contactinduced state change mechanism in the model accounts for spreading processes where the interaction with a given state can stimulate a node to adopt a state other than the one that is being interacted with. This extends the type of processes that can be modeled beyond epidemic spreading, and so the model also generalizes, for example, the Maki-Thompson rumor spreading Figure 13. A multiplex network. Different link types correspond to different layers in the multiplex network. We assume that transmission probabilities depend on the link type, i.e. each contagion or state propagates differently over each layer. This is depicted by coloring the contactinduced transition mechanism links differently for each separate layer. In the adaptation of the SI 1 I 2 S model for multiplex networks both contagions spread only on their respective layers. Hence, we have b