Skip to main content

Advertisement

Log in

Spreading dynamics on complex networks: a general stochastic approach

  • Published:
Journal of Mathematical Biology Aims and scope Submit manuscript

Abstract

Dynamics on networks is considered from the perspective of Markov stochastic processes. We partially describe the state of the system through network motifs and infer any missing data using the available information. This versatile approach is especially well adapted for modelling spreading processes and/or population dynamics. In particular, the generality of our framework and the fact that its assumptions are explicitly stated suggests that it could be used as a common ground for comparing existing epidemics models too complex for direct comparison, such as agent-based computer simulations. We provide many examples for the special cases of susceptible-infectious-susceptible and susceptible-infectious-removed dynamics (e.g., epidemics propagation) and we observe multiple situations where accurate results may be obtained at low computational cost. Our perspective reveals a subtle balance between the complex requirements of a realistic model and its basic assumptions.

This is a preview of subscription content, log in via an institution to check access.

Access this article

Price excludes VAT (USA)
Tax calculation will be finalised during checkout.

Instant access to the full article PDF.

Institutional subscriptions

Fig. 1
Fig. 2
Fig. 3
Fig. 4
Fig. 5
Fig. 6
Fig. 7
Fig. 8
Fig. 9

Similar content being viewed by others

Notes

  1. Computer science terminology is used for the sake of specificity. However, the reader should keep in mind that these simulations may be hypothetical: our goal is to replace them by a Level III model.

  2. While \(V\) should technically be deterministic, we may also perceive it as stochastic by considering pseudo-random number generators as “actually” random. In any case, these subtleties are of no concerns for our purpose.

  3. Remember that what we are ultimately interested in are the Level I open questions, so \(X(t)\) and \(Z(t)\) should be compared with respect to these concerns.

  4. The \(r_i^j\) may be positive, negative, or zero. For a given \(\mathbf{x}(t)\), if \(\mathbf{x}(t) \pm \mathbf{r}^j\) contains any negative entries, then we should have \(q_{j}^{\pm }\!\bigl ({\mathbf{x}(t)},{Y}\bigr ) = 0\) [i.e., \(\mathbf{x}(t+{{\mathrm{d}}}t)\) has zero probability to be negative].

  5. Notice that, although \(q_{j}^+\!\bigl ({\mathbf{x}(t)},{Y}\bigr )\) and \(q_{j}^-\!\bigl ({\mathbf{x}(t)},{Y}\bigr )\) tend to “cancel out” in the deterministic contribution [negative sign in (4a)], the stochastic contributions of transitions happening both forward and backwards instead “accumulate” [positive sign in (4d)].

  6. This is thus a special case of the configuration model, defined later, for which \(n_\kappa = N\) and \(n_{\kappa '} = 0 \ \forall \kappa ' \ne \kappa \).

  7. Compare (5b) with (5a): if a forward transition of type \(j\) occurs at time \(t\), we get \(\mathbf{x}(t+{{\mathrm{d}}}t) = \mathbf{x}(t) + \mathbf{r}^j\), so \(x_I(t+{{\mathrm{d}}}t) = x_I(t)+1\) and \(x_{S-I}(t+{{\mathrm{d}}}t) = x_{S-I}(t) + \kappa - 2j\). This event has probability \(q_{j}^+\!\bigl ({\mathbf{x}(t)},{Y}\bigr ){{\mathrm{d}}}t\) to occur during the time interval \([t,t+{{\mathrm{d}}}t)\).

  8. Recall that a backward transition of type \(j\) occurring at time \(t\) has the effect \(\mathbf{x}(t+{{\mathrm{d}}}t) = \mathbf{x}(t) - \mathbf{r}^j\). This event has probability \(q_{j}^-\!\bigl ({\mathbf{x}(t)},{Y}\bigr ){{\mathrm{d}}}t\) to occur during the time interval \([t,t+{{\mathrm{d}}}t)\).

  9. The actual distribution should be obtained by sampling without replacement, but for large networks we may make the approximation that the sampling is done with replacement, hence resulting in a binomial distribution.

  10. For example, a node of high degree is much more likely to acquire the infection, and is much more dangerous once it has acquired it, than a low degree node.

References

Download references

Acknowledgments

The research team acknowledges to the Canadian Institutes of Health Research (CIHR), the Natural Sciences and Engineering Research Council of Canada (NSERC) and the Fonds de recherche du Québec—Nature et technologies (FRQ–NT) for financial support. We are grateful to the anonymous referees that have helped improve our presentation.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Louis J. Dubé.

Appendices

Appendix A: More on motifs

We introduce new motifs and generalize those presented in the main text for intrinsic link types and directed (or semi-directed) networks; Table 1 summarizes the total number of motifs in each of these classes. We define \({\mathcal {K}}\) as the highest accessible node degree and we set \({\mathcal {D}} = 1\) for undirected, \({\mathcal {D}} = 2\) for directed and \({\mathcal {D}}= 3\) for semi-directed networks.

Table 1 Number of motifs in selected classes

Of high practical importance is the fact that the entries of this table differ widely in their scaling behaviours. For example, the on-the-fly model (Marceau et al. 2011) for two interacting SIR dynamics each propagating on their own network structure uses \({\mathcal {N}} = 9\) (Cartesian product of the sub-states of two SIR dynamics), \({\mathcal {L}} = 3\) (links of the first network alone, links of the second network alone and overlapping links) and \({\mathcal {D}} = 1\) (undirected network). Looking up in Table 1, we see that this requires of the order of \({\mathcal {K}}^3\) on-the-fly (degreed node) motifs, where \({\mathcal {K}}\) denotes the highest accessible node degree. By opposition, implementing a first neighbourhood version of this model would require of the order of \({\mathcal {K}}^{27}\) motifs!

1.1 A.1 Pair motifs

We note \(\nu \mathop {-}\limits ^{\ell }\nu '\) (resp. \(\nu \mathop {\rightarrow }\limits ^{\ell }\nu '\)) an undirected (resp. directed) pair motif formed of a state \(\nu \) node linked through a state \(\ell \) link to a state \(\nu '\) node. In the case of directed motifs, the direction of the arrow usually specifies the “strongest causal effect” of this asymmetric interaction (although this needs not be the case). All the undirected and directed motifs are possible in a semi-directed network. We may omit the index \(\ell \) over the links when \({\mathcal {L}} = 1\).

1.2 A.2 Triple motifs

We note \(\nu \mathop {-}\limits ^{\ell }\nu '\mathop {-}\limits ^{\ell '}\nu ''\) an undirected triple motif formed of a state \(\nu \) node linked through a state \(\ell \) link to a state \(\nu '\) node itself linked through a state \(\ell '\) link to a state \(\nu ''\) node. As for pair motifs, each of these nodes may have other neighbours than those that are explicitly specified. The notation directly generalizes to directed (e.g., \(\nu \mathop {\rightarrow }\limits ^{\ell }\nu '\mathop {\leftarrow }\limits ^{\ell '}\nu ''\)) and semi-directed (e.g., \(\nu \mathop {-}\limits ^{\ell }\nu '\mathop {\rightarrow }\limits ^{\ell '}\nu ''\)) triple motifs.

The term “\(2\)-star” is often used to refer to a triple motif for which both extremities (e.g., the nodes of state \(\nu \) and \(\nu ''\) in the motif \(\nu \mathop {-}\limits ^{\ell }\nu '\mathop {-}\limits ^{\ell '}\nu ''\)) are explicitly forbidden to be neighbours. In models that also use triangle motifs, \(2\)-star motifs may explicit the absence of the last link that would form a triangle. Another common use of triple motifs comes in the inference process of (usually deterministic) pair based models.

1.3 A.3 Triangle motifs and other small subnetworks

Three nodes that are all neighbours of each other form a triangle motif. An horizontal bracket represents the additional link that would be missing in the analogous triple motif, e.g.,

for an undirected network.

Triangle motifs are usually considered in models that should account for clustering. Their number may either directly be tracked in the state vector \(\mathbf{x}(t)\) (House et al. 2009) or their implicit presence (stated in \(Y\), e.g., through a clustering coefficient) may be accounted for in the inference process (Keeling et al. 1997).

The same notation may be generalized to other motifs consisting of small subnetworks, e.g., square motifs (House et al. 2009).

1.4 A.4 Clique motifs

A vague definition of a clique motif is “a subgroup of nodes that share more links among themselves than what could be expected otherwise for the same number of randomly selected nodes”. In applications, one may refine this definition according to the specificities of the problem at hand, e.g., “an Erdős-Rényi subnetwork (link probability \(p\)) of \(n_S\) susceptible nodes and \(n_I\) infectious nodes”. Clique motifs are usually considered in models that should account for community structure (Hébert-Dufresne et al. 2010).

1.5 A.5 Degreed motifs

A degreed motif is a motif for which we know the degree of all the nodes forming the motif: a degreed node motif is a node of specified state and degree; a degreed link motif is two nodes of specified state and degree that are known to be neighbours; etc. The in-degree and out-degree are both specified in directed and semi-directed networks; the latter cases also specify undirected degrees. Likewise, degrees pertaining to different types of links are specified independently.

In the same way that pair motifs are usually combined with node motifs, degreed pair motifs are usually combined with degreed node motifs (Eames and Keeling 2002). Note that the on-the-fly motifs \(\nu \Lambda _{1}({\kappa })\) presented Sect. 3.4 can be understood as a special case of degreed node motifs where the degree is replaced by the “degree to unknown nodes”.

1.6 A.6 \(n\)th neighbourhood motifs

Similarly to degreed motifs, an \(n\)th neighbourhood motif is a motif for which we specify the state of all the \(n\)th neighbours of the nodes forming the original motif. Hence, the notation \(\underline{\nu } \varGamma _{1}({\mathbf{k}})\) \(\bigl [\)resp. \(\underline{\nu } \varGamma _2({\mathbf{K}})\) \(\bigr ]\)of the main text corresponds to a first (resp. second) neighbourhood node motif. These concepts are directly generalizable to types of links and to directed or semi-directed networks.

First neighbourhood node motifs can be understood as tracking the correlation between the state of a node and the state of all its neighbours. By opposition, degreed pair motifs track the correlation between the state and degree of two neighbouring nodes. While similar information may be obtained from both motif classes, a model based on one may perform better than a model based on the other depending on the characteristics of the full system to be modelled.

1.7 A.7 Coarse-grained degree and/or neighbourhood

Not all entries of Table 1 depend on the maximal degree \({\mathcal {K}}\), but those that do quickly increase for large \({\mathcal {K}}\). This is problematic since many real-world systems contain high degree nodes.

However, one may overcome this limitation by coarse-graining degrees through ranges: the range containing the degree of a node is specified instead of the degree itself. For example, given the ranges

$$\begin{aligned} \underbrace{[0,0]}_{\text {range 0}}, \underbrace{[1,1]}_{\text {range 1}}, \underbrace{[2,3]}_{\text {range 2}}, \underbrace{[4,7]}_{\text {range 3}}, \underbrace{[8,15]}_{\text {range 4}}, \underbrace{[16,31]}_{\text {range 5}} \text { and } \underbrace{[32,63]}_{\text {range 6}}, \end{aligned}$$

we would say of a degree \(23\) node that its degree lies within range \(5\). Hence, for the purpose of evaluating the number of motifs in Table 1, one should here use \({\mathcal {K}} = 6\) (one less than the total number of ranges) instead of \({\mathcal {K}} = 63\) (highest representable degree).

While the previous example used powers of \(2\) for simplicity, a slower increase is probably desirable in most applications. However, since the number of neighbourhood and degreed motifs strongly depends on \({\mathcal {K}}\), even the slightest reduction in this number may be significant. Note that this coarse-graining method is of particular interest when the real-world data used to calibrate the model is already coarse-grained, which is commonly the case for census data.

Appendix B: Deterministic solution of on-the-fly SIR model

This section provides the deterministic solution of (9). We first rewrite (4b) for the specific case of (9)

$$\begin{aligned} \frac{{{\mathrm{d}}}\varvec{\mu }(t)}{{{\mathrm{d}}}t} = \sum _{\nu } \sum _{\kappa } \sum _{\kappa '} \mathbf{r}^{{\mathcal {I}}\nu \kappa \kappa '} q_{{\mathcal {I}}\nu \kappa \kappa '}^+\!\bigl ({\varvec{\mu }(t)},{Y}\bigr ) + \sum _{\kappa } \mathbf{r}^{{\mathcal {R}}\kappa } q_{{\mathcal {R}}\kappa }^+\!\bigl ({\varvec{\mu }(t)},{Y}\bigr ) \end{aligned}$$
(10)

then collect the contributions to \(\mu _{S\Lambda _{1}({\kappa })}\), \(\mu _{I\Lambda _{1}({\kappa })}\) and \(\mu _{R\Lambda _{1}({\kappa })}\) (dropping the functional dependencies for brevity)

$$\begin{aligned} \frac{{{\mathrm{d}}}\mu _{S\Lambda _{1}({\kappa })}}{{{\mathrm{d}}}t}&= -\sum _{\kappa '} q_{{\mathcal {I}}S\kappa \kappa '}^+ \end{aligned}$$
(11a)
$$\begin{aligned} \frac{{{\mathrm{d}}}\mu _{I\Lambda _{1}({\kappa })}}{{{\mathrm{d}}}t}&= \sum _{\kappa '} \left( q_{{\mathcal {I}}S(\kappa +1)\kappa '}^+ + q_{{\mathcal {I}}I(\kappa +1)\kappa '}^+ - q_{{\mathcal {I}}I\kappa \kappa '}^+ \right) \nonumber \\&+ \sum _{\nu } \sum _{\kappa '} \left( q_{{\mathcal {I}}\nu \kappa '(\kappa +1)}^+ - q_{{\mathcal {I}}I\kappa '\kappa }^+ \right) - q_{{\mathcal {R}}\kappa }^+ \end{aligned}$$
(11b)
$$\begin{aligned} \frac{{{\mathrm{d}}}\mu _{R\Lambda _{1}({\kappa })}}{{{\mathrm{d}}}t}&= \sum _{\kappa '} \left( q_{{\mathcal {I}}R(\kappa +1)\kappa '}^+ - q_{{\mathcal {I}}R\kappa \kappa '}^+ \right) + q_{{\mathcal {R}}\kappa }^+. \end{aligned}$$
(11c)

Using the definitions

$$\begin{aligned} \lambda = \sum _\kappa \kappa \mu _{I\Lambda _{1}({\kappa })} \quad \text {and} \quad \omega = \sum _\nu \sum _\kappa \kappa \mu _{\nu \Lambda _{1}({\kappa })} \end{aligned}$$
(12)

where \(\lambda \) is the total number of stubs belonging to infectious nodes and \(\omega \) is the total number of stubs in the system, (11) becomes

$$\begin{aligned} \frac{{{\mathrm{d}}}\mu _{S\Lambda _{1}({\kappa })}}{{{\mathrm{d}}}t}&= - \frac{\beta \lambda \kappa \mu _{S\Lambda _{1}({\kappa })}}{\omega } \end{aligned}$$
(13a)
$$\begin{aligned} \frac{{{\mathrm{d}}}\mu _{I\Lambda _{1}({\kappa })}}{{{\mathrm{d}}}t}&= \frac{\beta \lambda (\kappa + 1) \mu _{S\Lambda _{1}({\kappa + 1})}}{\omega } - \alpha \mu _{I\Lambda _{1}({\kappa })}\nonumber \\&+ \beta \left( 1 + \frac{\lambda }{\omega } \right) \Bigl ( (\kappa + 1) \mu _{I\Lambda _{1}({\kappa + 1})} - \kappa \mu _{I\Lambda _{1}({\kappa })} \Bigr ) \end{aligned}$$
(13b)
$$\begin{aligned} \frac{{{\mathrm{d}}}\mu _{R\Lambda _{1}({\kappa })}}{{{\mathrm{d}}}t}&= \frac{\beta \lambda }{\omega } \Bigl ( (\kappa + 1) \mu _{R\Lambda _{1}({\kappa + 1})} - \kappa \mu _{R\Lambda _{1}({\kappa })} \Bigr ) + \alpha \mu _{I\Lambda _{1}({\kappa })}. \end{aligned}$$
(13c)

Note that (9f) has been approximated by dropping the Kronecker delta in the numerator and the \(-1\) in the denominator.

We now consider the evolution of the total number of stubs \(\omega \) by summing the contributions from (13)

$$\begin{aligned} \frac{{{\mathrm{d}}}\omega }{{{\mathrm{d}}}t} = \sum _\nu \sum _\kappa \kappa \frac{{{\mathrm{d}}}\mu _{\nu \Lambda _{1}({\kappa })}}{{{\mathrm{d}}}t} = -2\beta \lambda . \end{aligned}$$
(14)

One may understand (14) as “during the time interval \([t,t+{{\mathrm{d}}}t)\), each one of the \(\lambda \) stubs belonging to infectious nodes have probability \(\beta {{\mathrm{d}}}t\) to be paired to another stub, thus causing a decrease by \(2\) of \(\omega \). Noting \(\omega _0 = \omega (0)\) the total number of stubs in the initial condition, we introduce the change of variable

$$\begin{aligned} \theta = \sqrt{ \frac{\omega }{\omega _0} } \quad \text {such that} \quad \frac{{{\mathrm{d}}}\theta }{{{\mathrm{d}}}t} = -\frac{\beta \lambda }{\theta \omega _0} \end{aligned}$$
(15)

Notice that \(t = 0\) corresponds to \(\theta = 1\) and that \(\theta \) decreases with time. Using this change of variable in (13a) gives

$$\begin{aligned} \frac{{{\mathrm{d}}}\mu _{S\Lambda _{1}({\kappa })}}{{{\mathrm{d}}}\theta } = \frac{\kappa \mu _{S\Lambda _{1}({\kappa })}}{\theta } \end{aligned}$$
(16)

which has the solution

$$\begin{aligned} \mu _{S\Lambda _{1}({\kappa })}(t) = x_{S\Lambda _{1}({\kappa })}(0) \bigl (\theta (t)\bigr )^\kappa \end{aligned}$$
(17)

using the initial condition \(\mu _{S\Lambda _{1}({\kappa })}(0) = x_{S\Lambda _{1}({\kappa })}(0)\).

For convenience, we define

$$\begin{aligned} f(\theta ) = \sum _{\kappa } x_{S\Lambda _{1}({\kappa })}(0) \theta ^\kappa . \end{aligned}$$
(18)

Noticing that

$$\begin{aligned} \sum _{\kappa } \kappa (\kappa + 1) \mu _{S\Lambda _{1}({\kappa +1})}&= \theta ^2 \sum _{\kappa } (\kappa - 1) \kappa x_{S\Lambda _{1}({\kappa })}(0) \theta ^{\kappa -2} = \theta ^2 f''(\theta ) \quad , \end{aligned}$$
(19)

we obtain the evolution of the total number \(\lambda \) of stubs belonging to infectious nodes by summing the contributions (13b)

$$\begin{aligned} \frac{{{\mathrm{d}}}\lambda }{{{\mathrm{d}}}t}&= \sum _\kappa \kappa \frac{{{\mathrm{d}}}\mu _{I\Lambda _{1}({\kappa })}}{{{\mathrm{d}}}t} = \frac{\beta \lambda \theta ^2 f''(\theta )}{\omega } - \beta \lambda \left( 1 + \frac{\lambda }{\omega } \right) - \alpha \lambda . \end{aligned}$$
(20)

Again using the change of variable (15), we get

$$\begin{aligned} \frac{{{\mathrm{d}}}\lambda }{{{\mathrm{d}}}\theta } = \frac{\lambda }{\theta } + \theta \omega _0 \left( 1 + \frac{\alpha }{\beta } \right) - \theta f''(\theta ) \end{aligned}$$
(21)

which has the solution

$$\begin{aligned} \lambda = \theta ^2 \omega _0 \left( 1 + \frac{\alpha }{\beta } \right) - \theta \omega _0 \frac{\alpha }{\beta } - \theta f'(\theta ) \end{aligned}$$
(22)

for an initial condition without removed nodes \(\bigl [\) i.e., \(\lambda (0) = \omega _0 - f'(1)\) \(\bigr ]\).

Using this solution in (15) gives

$$\begin{aligned} \frac{{{\mathrm{d}}}\theta }{{{\mathrm{d}}}t} = - \beta \theta + \alpha \left( 1 - \theta \right) + \beta \frac{f'(\theta )}{\omega _0} \end{aligned}$$
(23)

whose solution provides \(\theta (t)\). Using

$$\begin{aligned} S(t)&= f\bigl (\theta (t)\bigr ) \end{aligned}$$
(24a)
$$\begin{aligned} I(t)&= N - S(t) - R(t) \end{aligned}$$
(24b)
$$\begin{aligned} \frac{{{\mathrm{d}}}R(t)}{{{\mathrm{d}}}t}&= \alpha I(t), \end{aligned}$$
(24c)

we finally obtain the total number \(S = \sum _\kappa \mu _{S\Lambda _{1}({\kappa })}\) of susceptibles, \(I = \sum _\kappa \mu _{I\Lambda _{1}({\kappa })}\) of infectious and \(R = \sum _\kappa \mu _{R\Lambda _{1}({\kappa })}\) of removed nodes at any given time \(t\). A direct application of (17)–(18) provides (24a), conservation of the nodes provides (24b) and using the definitions of \(I(t)\) and \(R(t)\) in (13c) provides (24c). Although obtained differently, this solution corresponds to that of the pair-based SIR model presented in Volz (2008) and Miller (2010).

Rights and permissions

Reprints and permissions

About this article

Cite this article

Noël, PA., Allard, A., Hébert-Dufresne, L. et al. Spreading dynamics on complex networks: a general stochastic approach. J. Math. Biol. 69, 1627–1660 (2014). https://doi.org/10.1007/s00285-013-0744-9

Download citation

  • Received:

  • Published:

  • Issue Date:

  • DOI: https://doi.org/10.1007/s00285-013-0744-9

Keywords

Mathematics Subject Classification (2000)

Navigation