Lumping of Degree-Based Mean Field and Pair Approximation Equations for Multi-State Contact Processes

TL;DR

This paper introduces a method to simplify and accelerate the analysis of multi-state contact processes on large networks by clustering nodes with similar degrees, significantly reducing computational complexity while maintaining accuracy.

Contribution

It extends degree-based mean field and pair approximation methods to multi-state processes and provides an automatic clustering approach for efficient analysis.

Findings

High compression of equations achieved

Significant reduction in computational time

Minimal loss of accuracy in results

Abstract

Contact processes form a large and highly interesting class of dynamic processes on networks, including epidemic and information spreading. While devising stochastic models of such processes is relatively easy, analyzing them is very challenging from a computational point of view, particularly for large networks appearing in real applications. One strategy to reduce the complexity of their analysis is to rely on approximations, often in terms of a set of differential equations capturing the evolution of a random node, distinguishing nodes with different topological contexts (i.e., different degrees of different neighborhoods), like degree-based mean field (DBMF), approximate master equation (AME), or pair approximation (PA). The number of differential equations so obtained is typically proportional to the maximum degree kmax of the network, which is much smaller than the size of the master equation of the underlying stochastic model, yet numerically solving these equations can still be problematic for large kmax. In this paper, we extend AME and PA, which has been proposed only for the binary state case, to a multi-state setting and provide an aggregation procedure that clusters together nodes having similar degrees, treating those in the same cluster as indistinguishable, thus reducing the number of equations while preserving an accurate description of global observables of interest. We also provide an automatic way to build such equations and to identify a small number of degree clusters that give accurate results. The method is tested on several case studies, where it shows a high level of compression and a reduction of computational time of several orders of magnitude for large networks, with minimal loss in accuracy.

Tables

Table 1. Table 1: SIR model: Running times (in seconds) of the full equations ( T full subscript 𝑇 full T_{\mathrm{full}} ) vs. running times of the lumped equations ( T lumped subscript 𝑇 lumped T_{\mathrm{lumped}} ) for all network sizes.

$k_{\max }$	DBMF($T_{\mathrm{full}}$)	DBMF($T_{\mathrm{lumped}}$)	PA($T_{\mathrm{full}}$)	PA($T_{\mathrm{lumped}}$)
${10}^2$	3	2	1	0.07
${10}^3$	5	3	9	0.06
${10}^4$	10	2	126	0.32
${10}^5$	228	3	4428	3
${10}^6$	18177	3	449280	12

Table 2. Table 2: SIR model: Number of bins returned by Alg.  1 .

$k_{\max }$	DBMF	DBMF $ϵ$	PA	PA $ϵ$
${10}^2$	20	$3.80\cdot {10}^{-7}$	20	$2.70\cdot {10}^{-8}$
${10}^3$	20	$6.05\cdot {10}^{-5}$	20	$7.95\cdot {10}^{-6}$
${10}^4$	55	$1.04\cdot {10}^{-4}$	56	$1.80\cdot {10}^{-4}$
${10}^5$	37	$4.47\cdot {10}^{-4}$	37	$4.52\cdot {10}^{-4}$
${10}^6$	35	$5.87\cdot {10}^{-4}$	35	$6.18\cdot {10}^{-4}$

Equations

$$\frac{\partial x_s}{\partial t}$$

$$\frac{\partial x_{s,k}}{\partial t}$$

$$p_k[s^{\prime }]=p[s^{\prime }]=\sum _{k^{\prime }\in K}\frac{P(k^{\prime })k^{\prime }{x}_{s^{\prime },k^{\prime }}}{⟨k⟩},$$

$$\frac{\partial x_{s,k,m}}{\partial t}=\sum _{R_{I,i}\in R_I^{s^{+}}}{\lambda }_i{x}_{r_i,k,m}-\sum _{R_{I,i}\in R_I^{s^{-}}}{\lambda }_i{x}_{s,k,m}$$

$${\beta }^{s{s}^{\prime }\to s{s}^{\prime \prime }}=\frac{1}{⟨{\sum }_{m\in {\mathbb{M}}_k}m[s]x_{s^{\prime },k,m}⟩}\cdot$$

$$\displaystyle\left\langle\sum_{\mathbb{m}\in\mathbb{M}_{k}}\!\!\mathbb{m}[s]x_{s^{\prime},k,\mathbb{m}}\Big{(}\sum_{R_{I,i}\in R^{s^{\prime}\rightarrow s^{\prime\prime}}_{I}}\!\!\!\!\lambda_{i}+\!\!\!\sum_{R_{C,j}\in R^{s^{\prime}\rightarrow s^{\prime\prime}}_{C}}\!\!\!\!\!\lambda_{j}\mathbb{m}[r_{j_{2}}]\Big{)}\right\rangle,$$

$$\frac{\partial x_{s,k}}{\partial t}=\sum _{R_{I,i}\in R_I^{s^{+}}}{\lambda }_i{x}_{r_i,k}-\sum _{R_{I,i}\in R_I^{s^{-}}}{\lambda }_i{x}_{s,k}$$

$$+\sum _{R_{C,j}\in R_C^{s^{+}}}{\lambda }_j{x}_{r_{j_1},k}{\mu }_{r_{j_1},k}(r_{j_2})-\sum _{R_{C,j}\in R_C^{s^{-}}}{\lambda }_j{x}_{s,k}{\mu }_{s,k}(r_{j_2})$$

$$\frac{\partial p_{s,k}[s_n]}{\partial t}=-\frac{\partial x_{s,k}}{\partial t}\frac{p_{s,k}[s_n]}{x_{s,k}}$$

$$+\sum _{R_{I,i}\in R_I^{s^{+}}}\frac{{\lambda }_i{x}_{r_i,k}}{x_{s,k}k}{\mu }_{r_{i,k}}(s_n)-\sum _{R_{I,i}\in R_I^{s^{-}}}\frac{{\lambda }_i}{k}{\mu }_{s,k}(s_n)$$

$$+\sum _{R_{C,j}\in R_C^{s^{+}}}\frac{{\lambda }_j{x}_{r_{j_1},k}}{x_{s,k}k}{\mu }_{r_{j_1},k}(r_{j_2},s_n)-\sum _{R_{C,j}\in R_C^{s^{-}}}\frac{{\lambda }_j}{k}{\mu }_{s,k}(r_{j_2},s_n)$$

$$+\sum _{\begin{array}{c}{s}^{\prime }\in S\\ s^{\prime }\ne s_n\end{array}}{\beta }^{s{s}^{\prime }\to s{s}_n}{p}_{s,k}[s^{\prime }]-\sum _{\begin{array}{c}{s}^{\prime }\in S\\ s^{\prime }\ne s_n\end{array}}{\beta }^{s{s}_n\to s{s}^{\prime }}{p}_{s,k}[s_n],$$

$$\displaystyle\beta^{ss^{\prime}\rightarrow ss^{\prime\prime}}=1/{\left\langle x_{s^{\prime},k}{\mu}_{s^{\prime},k}(s)\right\rangle}\cdot\big{(}\big{\langle}x_{s^{\prime},k}\sum_{R^{s^{\prime}\rightarrow s^{\prime\prime}}_{I}}\lambda_{i}{\mu}_{s^{\prime},k}(s)\big{\rangle}$$

$$\displaystyle+\ \big{\langle}x_{s^{\prime},k}\sum_{R^{s^{\prime}\rightarrow s^{\prime\prime}}_{C}}\lambda_{j}{\mu}_{s^{\prime},k}(s,r_{j_{2}})\big{\rangle}\big{)}.$$

$$\frac{\partial }{\partial t}\sum _{k\in b}{x}_{s,k}P(k|b)=$$

$$\sum _{R_{I,i}\in R_I^{s^{+}}}{\lambda }_i\sum _{k\in b}{x}_{r_i,k}P(k|b)-\sum _{R_{I,i}\in R_I^{s^{-}}}{\lambda }_i\sum _{k\in b}{x}_{s,k}P(k|b)$$

$$+\sum _{R_{C,j}\in R_C^{s^{+}}}{\lambda }_{j}p[r_{j_2}]\sum _{k\in b}{x}_{r_{j_1,k}}kP(k|b)$$

$$-\sum _{R_{C,j}\in R_C^{s^{-}}}{\lambda }_{j}p[r_{j_2}]\sum _{k\in b}{x}_{s,k}kP(k|b)$$

$$k\in b$$

$$\frac{\partial }{\partial t}{x}_{s,b}=\sum _{R_{I,i}\in R_I^{s^{+}}}{\lambda }_i{x}_{r_i,b}-\sum _{R_{I,i}\in R_I^{s^{-}}}{\lambda }_i{x}_{s,b}$$

$$\begin{array}{rl}p[s_n]& =\frac{1}{⟨k⟩}\sum _{k^{\prime }\in K}P(k^{\prime })k^{\prime }{x}_{s_n,k^{\prime }}\\ & =\frac{1}{⟨k⟩}\sum _{b\in B}\sum _{k^{\prime }\in b}P(k^{\prime })k^{\prime }{x}_{s_n,k^{\prime }}\\ & =\frac{1}{⟨k⟩}\sum _{b\in B}{x}_{s_n,b}P(b)\sum _{k^{\prime }\in b}P(k^{\prime }|b)k^{\prime }\\ & =\frac{1}{⟨k⟩}\sum _{b\in B}{x}_{s_n,b}P(b){⟨k⟩}_b.\end{array}$$

$$\frac{\partial }{\partial t}\sum _{k\in b}{x}_{s,k}P(k|b)=$$

$$\sum _{R_{I,i}\in R_I^{s^{+}}}{\lambda }_i\sum _{k\in b}{x}_{r_i,k}P(k|b)-\sum _{R_{I,i}\in R_I^{s^{-}}}{\lambda }_i\sum _{k\in b}{x}_{s,k}P(k|b)$$

$$+\sum _{R_{C,j}\in R_C^{s^{+}}}{\lambda }_j\sum _{k\in b}{x}_{r_{j_1},k}{\mu }_{r_{j_1},k}(r_{j_2})P(k|b)$$

$$-\sum _{R_{C,j}\in R_C^{s^{-}}}{\lambda }_j\sum _{k\in b}{x}_{s,k}{\mu }_{s,k}(r_{j_2})P(k|b),$$

$${\mu }_{s,b}(s_n)$$

$$\frac{\partial }{\partial t}{x}_{s,b}$$

$$\frac{\partial }{\partial t}\sum _{k\in b}{p}_{s,k}[s_n]P(k|b)=-\sum _{k\in b}\frac{\partial }{\partial t}{x}_{s,k}\frac{p_{s,k}[s_n]}{x_{s,k}}P(k|b)$$

$$\frac{\partial }{\partial t}{p}_{s,b}[s_n]=-\frac{\partial }{\partial t}{x}_{s,b}\frac{p_{s,b}[s_n]}{x_{s,b}}+\sum _{R_{I,i}\in R_I^{s^{+}}}{\lambda }_i\frac{x_{r_i,b}}{x_{s,b}}{p}_{r{}_i,b}[s_n]$$

$$\sum _{k\in b}\frac{1}{k}{\mu }_{s,k}(r_{j_2},s_n)P(k|b)={⟨k-1⟩}_b{p}_{s,b}[r_{j_2}]p_{s,b}[s_n]$$

$$\sum _{k\in b}\frac{1}{k}{\mu }_{s,k}(r_{j_2},s_n)P(k|b)=p_{s,b}[s_n]+p_{s,b}[s_n{]}^2{⟨k-1⟩}_b.$$

$${\mu }_{s,b}(s_n,s_n^{\prime })=\sum _{k\in b}{\mu }_{s,k}(s_n,s_n^{\prime })P(k|b),$$

$${\beta }_B^{s{s}^{\prime }\to s{s}^{\prime \prime }}=1/{⟨x_{s^{\prime },b}{\mu }_{s^{\prime },b}(s)⟩}_b$$

$$\displaystyle\big{(}{\big{\langle}x_{s^{\prime},b}\sum_{R^{s^{\prime}\rightarrow s^{\prime\prime}}_{I}}\lambda_{i}{\mu}_{s^{\prime},b}(s)\big{\rangle}_{b}+\big{\langle}x_{s^{\prime},b}\sum_{R^{s^{\prime}\rightarrow s^{\prime\prime}}_{C}}\lambda_{j}{\mu}_{s^{\prime},b}(s,r_{j_{2}})\big{\rangle}_{b}}\big{)}.$$

$$\begin{array}[]{rcl}\epsilon_{B}(t)&=&||x_{s}(t)-\bar{x}_{s}(t)||\\[4.30554pt] &=&||\sum_{b\,\in B}\sum_{k\in b}x_{s,k}(t)P(k)-x_{s,b}(t)P(k)||\\[4.30554pt] &=&||\sum_{b\,\in B}\sum_{k\in b}P(k)\underbrace{(x_{s,k}(t)-x_{s,b}(t))}_{\epsilon_{k,b}(t)}||.\end{array}$$

$$\displaystyle\begin{split}\text{dist}(b_{i},b_{j})&=P(b_{i})+P(b_{j})+\alpha\;\text{dist}_{hmg}(b_{i},b_{j}),\\ \text{dist}_{hmg}(b_{i},b_{j})&=\frac{\big{|}\bar{k}_{b_{i}}-\bar{k}_{b_{j}}\big{|}}{\big{(}\bar{k}_{b_{i}}+\bar{k}_{b_{j}}\big{)}/2},\end{split}$$

Full text

††preprint: APS/123-QED

I Introduction

Dynamic processes on complex networks have attracted a lot of interest in recent years Albert and Barabási (2002); Newman (2003); Strogatz (2001); Palla et al. (2005); Boccaletti et al. (2006); Barrat et al. (2008); Dorogovtsev et al. (2008); Pastor-Satorras et al. (2015). Their behavior is complex and their dynamics is strongly influenced by the network topology through interactions between connected nodes. Even more challenging is the study of adaptive networks, where both the node state and the network topology coevolve in a feedback relationship Gross et al. (2006); Gross and Blasius (2008); Marceau et al. (2010); Van Segbroeck et al. (2010); Wu et al. (2012); Wang et al. (2015).

Among dynamic processes happening on networks, a prominent role is played by contact processes, which include information spreading and epidemic models Barrat et al. (2008); Mata et al. (2014); Gómez et al. (2011); Juhász et al. (2012); Boguñá et al. (2009); Lee et al. (2013); Pastor-Satorras et al. (2015). In a contact process, each node of the network can be in one of several states, and neighboring nodes can interact and change their state. Interactions typically happen at random times, governed by an exponential distribution. The underlying stochastic process is therefore a discrete-state Markov jump process.

The analysis of such models is very challenging, because of the huge state space of the underlying process: each possible network configuration (i.e., each assignment of states to nodes) represents a different discrete (global) state of the process. The exponential blowup of the state space rules out the use of any numerical technique to compute probabilities, such as uniformization or finite projection methods Didier et al. (2009); Wolf et al. (2010), but also renders simulation quite expensive, particularly for very large networks. An alternative and viable approach for large networks is provided by mean-field-like approximations Bortolussi et al. (2013). There are several forms and variations of such approximations Pastor-Satorras and Vespignani (2001); Gleeson (2011); Cator and Van Mieghem (2012); Mata and Ferreira (2013); Demirel et al. (2014); Pastor-Satorras et al. (2015) for network models, all sharing the same basic idea, namely to treat all nodes with a similar local structure as indistinguishable, and write a differential equation for each class of nodes, modeling the evolution of a random node in time. Depending on how much information about the local structure is taken into account, we obtain more refined approximations, at the price of more equations to be solved. In the simplest case, we treat all nodes with the same local state as indistinguishable (mean-field (MF) approximation Bortolussi et al. (2013); Barrat et al. (2008)). Improvements can be obtained by grouping nodes with respect to their state and degree, called degree-based (described also as heterogeneous sometimes) mean-field (DBMF) Pastor-Satorras and Vespignani (2001); Castellano and Pastor-Satorras (2010); Gleeson et al. (2012), while an even better approximation is achieved by also representing explicitly also the state of neighboring nodes, as introduced by Gleeson’s approximate master equation (AME) for the binary state case Gleeson (2011). By implicitly relying on a binomial assumption, these equations can be simplified to yield the pair approximation (PA) Gleeson (2011). Although for many network topologies AME and PA give significantly more accurate results Gleeson (2013), solving the AME and PA, even for the binary state case considered by Gleeson, is much more expensive than DBMF.

For the general multi-state case that we consider here, the solution of AME and PA is usually infeasible for many real networks, as we show that the number of equations is $\binom{k_{\mathrm{max}}+|S|}{|S|-1}(k_{\mathrm{max}}+1)$ for AME and $(|S|^{2}+|S|)k_{\mathrm{max}}$ for PA, where $k_{\mathrm{max}}$ is the maximum degree of the network and $|S|$ is the number of possible states of a node.

Typically, when analyzing dynamic processes on complex networks, the interest is in emergent behaviors that manifest at a collective level. Hence, aggregated quantities like the fraction of nodes in a certain state are typically monitored. However, the network topology, even if marginalized in such observables, highly influences their dynamic properties. Often, the necessary granularity of the analysis is on the level of DBMF or PA. Here, we show that the level of detail of these approaches can be further decreased. Our lumping approach shows that for an accurate approximation of the mean behavior of the nodes, the exact degree and state of the neighborhood need not be monitored. Instead, aggregated information is sufficient, which strongly lowers the complexity of the analysis. We decrease the number of equations by several orders of magnitude without noticable loss in precision compared to PA or DBMF. Moreover, we are able to set up an automatic procedure to construct such compressed equations, which avoids altogether the solution of the full system of equations. Computational gains are impressive, up to five orders of magnitude for large networks (large $k_{\mathrm{max}}$). This increased performance can be extremely useful when equations have to be solved many times, for instance for parameter sweeping or for parameter identification tasks.

In this paper, we consider AME and PA for the case of $|S|$ local states, providing explicitly the AME and the derivation of PA in a similar fashion as in Fennell (2015). For the special case of contact processes, we arrive at an equation that involves the first and second moments of the multinomial distribution, which renders the equations computationally more tractable compared to Fennell (2015).

Our major contributions are the following: (i) we derive a lumped set of equations for the multi-state DBMF and PA equations and propose an algorithm to effectively construct this set; and (ii) we provide a numerical study on several non-trivial examples, showing that our approach performs well for several classes of degree distributions, that the reduction becomes more and more effective as the maximum network degree increases, and that non-linear effects in the dynamics are captured well by the reduced equations. Furthermore, (iii) we produced and made publicly available an open-source Python implementation of our method, taking a compact description of the model as input and returning Python files containing the lumped (and the full) DBMF and PA equations as output, together with procedures to integrate them in time and to plot the results.

The paper is further organized as follows: In Section II, we introduce contact processes and several approximations, including a generalization of AME and PA for $|S|$ local states. We derive lumped equations, given a clustering of degrees of the network, in Section III. In Section IV we propose a clustering strategy that yields the lumped equations resulting in an accurate solution, and present several case studies. Finally, in Section V we draw final conclusions.

II Approximate Equations

In this section, we first introduce contact processes on networks, and then discuss several mean-field approaches to approximate their dynamics. These encompass mean-field, degree-based mean-field, approximate master equation, and pair approximation. The latter two methods are generalized from Gleeson (2011) to models with $|S|$ local states.

II.1 Contact processes on networks

A contact process among individuals on a finite graph $G$ can be described as an interaction model $M=\{S,R_{I},R_{C},G\},$ where $S=\{s_{1},s_{2},\ldots,s_{|S|}\}$ contains the possible states of an individual, $R_{I}=\{R_{I,1},\ldots,R_{I,N_{I}}\}$ is the set of independent rules and $R_{C}=\{R_{C,1},\ldots,R_{C,N_{C}}\}$ is the set of the contact rules among individuals.

An independent rule $R_{I,i}\in R_{I},\;i\in\{1,\ldots,N_{I}\}$ models a spontaneous state change of a node of the graph $G$ from $r_{i}$ to $p_{i}$, where both $r_{i},p_{i}\in S$ and $r_{i}\neq p_{i}$. One instance of such a state change happens at an average rate of $\lambda_{i}$ per node. Independent rules are illustrated in Figure 1.

A contact rule $R_{C,j}\in R_{C},\;j\in\{1,\ldots,N_{C}\}$ describes the interaction of two neighboring nodes, respectively in states $r_{j_{1}},r_{j_{2}}\in S$, and changes the state of one such node, say state $r_{j_{1}}$, to $p_{j_{1}}\in S$, while the state of the other node remains the same, i.e, $r_{j_{2}}=p_{j_{2}}$. The general form of a contact rule is depicted in Figure 2. By abuse of notation, we define the average rate of this event per interacting edge as $\lambda_{j}$, i.e., we omit subscripts $I$ and $C$ to distinguish rates of independent and contact rules.

To express different variants of approximate master equations for an individual in state $s\in S$, we consider the subset of either the independent or the contact rules that affect the individual’s state. Let us define $R^{s^{+}}_{I}=\{R_{I,i}\in R_{I}\mid p_{i}=s\}$ as the subset of independent rules according to which a node spontaneously changes its state to $s.$ Similarly, $R^{s^{-}}_{I}=\{R_{I,i}\in R_{I}\mid r_{i}=s\}$ is the subset of independent rules where a node in state $s$ changes to some other state. For contact rules that affect the population of individuals in state $s$, we define $R^{s^{+}}_{C}=\{R_{C,j}\in R_{C}\mid p_{j_{1}}=s\}$ for those contact rules that add a new individual in state $s$ and $R^{s^{-}}_{C}=\{R_{C,j}\in R_{C}\mid r_{j_{1}}=s\}$ for those rules that result in a loss of an $s$ node.

The above contact process can be interpreted as a continuous-time Markov chain $\{X(t),t\geq 0\}$ with state space $L$, where each global state $g\in L$ corresponds to a labeled version of the graph $G$. The nonzero transition rates between the different global states are given by the rate of either an independent or a contact rule, where we assume that a global state change is triggered by at most one rule. Since the size of $L$ is $|S|^{N},$ where $N$ is the number of nodes in $G$, the analysis of $X$ is notoriously difficult. Typically, Monte Carlo simulations are carried out or equations for the approximation of its mean behavior are derived, since standard numerical approaches are infeasible.

II.2 Mean-Field Equations

One of the simplest assumptions to derive approximations of $X$ is that the nodes of the network are independent and equivalent. This is known as the mean-field assumption and gives rise to the well-known mean-field (MF) equation Barrat et al. (2008).

[TABLE]

In (II.2), $x_{s}(t)$ is the fraction of nodes in state $s$ at time $t$, $\left\langle k\right\rangle$ is the average degree of the network and $\lambda_{i},\lambda_{j}$ are the rates associated with the rules $R_{I,i},R_{C,j}$, respectively. The first two sums express the cases where a node changes its state to or from $s$ due to an independent rule and the last two sums represent every contact rule $R_{C,j}$ that either produces or consumes an $s$ node. The MF assumption is known to approximate well the dynamics of contact processes on either constant or Erdos-Renyi random graphs, where the degree of each node does not deviate much from the average degree Barrat et al. (2008).

If the network topology is slightly more complicated, which is the case for most real-world networks, then the MF equations cannot capture the dynamics of the spreading process anymore. In such networks, the infection probability of a node depends on the node’s degree. Consequently, the approximation has to take into account the network’s degree distribution $P$. This is particularly important for power law degree distributions, which are known to have a profound impact on the dynamics of contact processes Barrat et al. (2008); Broder et al. (2000); Newman (2001); Clauset et al. (2009); Pastor-Satorras et al. (2015). An approximation that distinguishes between the different node degrees in the network, while averaging over all nodes that have the same degree and state, is the degree-based mean-field (DBMF) Pastor-Satorras and Vespignani (2001). Let $K=\{1,\ldots,k_{\rm{max}}\}$ be the set of node degrees of the network and let $x_{s,k}$ be the fraction of nodes of degree $k$ that are in state $s\in S.$ We assume that the network has no nodes of degree zero, as they are irrelevant for the contact dynamics. Then, for general contact processes, we have for all $s\in S$ and all $k\in K$ that

[TABLE]

Here, $p_{k}[s^{\prime}]=Pr(s^{\prime}\mid k)$ denotes the probability that a random neighbor of a degree $k$ node is in state $s^{\prime}\in S$. We also let $p_{k}$ be the vector with entries $p_{k}[s^{\prime}]$, $s^{\prime}\in S$. If the network is uncorrelated, this probability is by definition independent of $k$ and for all $k\in K$ it holds (see also Appendix A):

[TABLE]

where ${\left\langle k\right\rangle}$ is the average degree. Hence, for uncorrelated networks the system of ODEs in (II.2) can be solved for some given degree distribution $P$ and initial conditions. Having computed $x_{s,k}$ for all $k\in K,$ we can determine $x_{s}=\sum_{k\in K}x_{s,k}P(k),$ which is usually the fraction of interest.

II.3 Approximate Master Equation and Pair Approximation

Mean-field theories (MF and DBMF) are known to be accurate for well-connected networks. However, they often perform poorly on sparse networks Gleeson (2011). A more accurate approximation is the approximate master equation (AME), which has been proposed for binary state contact processes, such as the SIS model of epidemic spreading Gleeson (2011); Marceau et al. (2010). Here we generalize the AME to contact processes with an arbitrary state space $S$.

For each $(s,k)$ node, i.e., each node that is in state $s$ and has degree $k,$ we define a neighbor vector $\mathbb{m}$ such that the entry $\mathbb{m}[s_{n}]:0\leq\mathbb{m}[s_{n}]\leq k$ contains the number of edges that connect the $(s,k)$ node to a node in state $s_{n}$, that is, all entries are non-negative integers and the sum of all entries is $k$. Let $\mathbb{M}_{k}$ be the set of all such vectors. For $x_{s,k,\mathbb{m}},$ i.e., the fraction of nodes of degree $k$ which are in state $s$ and have neighbor vector $\mathbb{m}$, we get for all $s\in S$, $k\in K$, and $\mathbb{m}\in\mathbb{M}_{k},$

[TABLE]

The first four terms in the above equation are similar to those in Eq. (II.2) and express the cases where we either add or delete an $(s,k,\mathbb{m})$ node via an independent or a contact rule. The difference now is that we can write down the precise number of edges between two reacting neighboring nodes in states $r_{j_{1}},r_{j_{2}}$ instead of using their average number as in (II.2).

Each term in the last two sums represents the case where one of the $k$ neighbors of an $(s,k)$ node changes its state from $s^{\prime}$ to $s^{\prime\prime}$. This change results in either obtaining a new $(s,k,\mathbb{m})$ node (first sum) or losing such a node (second sum).

The vector $\mathbb{m}^{\{s^{\prime+},s^{\prime\prime-}\}}$ is such that all entries are equal to those of $\mathbb{m}$, apart from the $s^{\prime}$-th entry, which is equal to $\mathbb{m}[s^{\prime}]+1$, and the $s^{\prime\prime}$-th entry, which is equal to $\mathbb{m}[s^{\prime\prime}]-1.$ In addition, $\beta^{ss^{\prime}\rightarrow ss^{\prime\prime}}$ is the average rate at which $ss^{\prime}$ edges change to $ss^{\prime\prime}$ edges, where $s,s^{\prime},s^{\prime\prime}\in S.$ To derive an expression for $\beta^{ss^{\prime}\rightarrow ss^{\prime\prime}}$ we define the sets $R^{s^{\prime}\rightarrow s^{\prime\prime}}_{I}=\{R_{I,i}\in R_{I}\mid r_{i}=s^{\prime}\wedge p_{i}=s^{\prime\prime}\}$ as the set of independent rules and $R^{s^{\prime}\rightarrow s^{\prime\prime}}_{C}=\{R_{C,j}\in R_{C}\mid r_{j_{1}}=s^{\prime}\wedge p_{j_{1}}=s^{\prime\prime}\}$ as the set of contact rules that result in a node change from $s^{\prime}$ to $s^{\prime\prime}$. Now, we define the rate $\beta^{ss^{\prime}\rightarrow ss^{\prime\prime}}$, which is independent of $k$, as

[TABLE]

where the above average is taken over all degrees of the network, i.e., $\left\langle f(k)\right\rangle=\sum_{k=1}^{k_{\mathrm{max}}}P(k)f(k).$

The number of differential equations that need to be solved for the general case of the AME is $|S|\sum_{k=0}^{k_{\mathrm{max}}}\binom{k+|S|-1}{|S|-1}=\binom{k_{\mathrm{max}}+|S|}{|S|-1}(k_{\mathrm{max}}+1).$ The binomial sum arises from the number of ways that, for a fixed degree $k,$ one can distribute $k$ unlabeled edges to $|S|$ states (see also Appendix B). For the special case of $|S|=2,$ the number of equations is $(k_{\mathrm{max}}+1)(k_{\mathrm{max}}+2)$ Gleeson (2011).

A coarser approximation, the pair approximation Gleeson (2011), is possible if instead of considering one equation for each vector $\mathbb{m}$ we assume a multinomial distribution for the number of edges that a node of degree $k$ in state $s$ has to all other states. This way, one needs to consider a (time-dependent) vector $p_{s,k}$ whose $i-$th entry $p_{s,k}[s_{i}]$ is defined as the probability that a random neighbor of an $(s,k)$ node is in state $s_{i}\in S$. Here we consider the general pair approximation for contact processes with $|S|$ states, generalizing the binary state approach of Gleeson (2011) similarly to Fennell (2015). A detailed derivation is provided in Appendix C. Assuming that $x_{s,k}>0$, we have for all $s\in S$, $k\in K$, $s_{n}\in S$,

[TABLE]

and

[TABLE]

where ${\mu}_{s,k}(s_{n})=\mathbb{E}\left[\mathbb{m}[s_{n}]\right]$ is the expected number of $s_{n}$ neighbors of an $(s,k)$ node. We assume that the number of $s_{n}$ neighbors is binomially distributed, i.e., the distribution is $\mathrm{\textbf{B}}(k,p_{s,k}[s_{n}])$. Furthermore, ${\mu}_{s,k}(s_{n},s^{\prime}_{n})=\mathbb{E}\left[\mathbb{m}[s_{n}]\mathbb{m}[s^{\prime}_{n}]\right]$ is the second mixed moment of the corresponding multinomial distribution with parameter vector $p_{s,k}.$ In a similar fashion, exploiting the multinomial assumption we can further simplify (see Appendix C) $\beta^{ss^{\prime}\rightarrow ss^{\prime\prime}}$ to

[TABLE]

The number of equations is now reduced to $(|S|^{2}+|S|)k_{\mathrm{max}}$, where $|S|$ is typically much smaller than $k_{\mathrm{max}}.$ In addition, the use of the moments of the multinomial distribution results in an enormous computational gain compared to the enumeration of huge sums over vectors when solving the above ODE system. For a detailed derivation of the above equations and a proof of their correctness we refer the reader to Appendix C.

As in the pair approximation for the binary case Gleeson (2011), the right hand side of Eq. (II.3) is undefined if $x_{s,k}=0$. However, it is easy to check that if $x_{s,k}\neq 0,$ then the system can never reach such a value, as the negative terms in Eq. (6) depend linearly on $x_{s,k}$. In the following, without loss of generality, we assume that $x_{s,k}\neq 0$ for all $s$ and $k$ at all times. An alternative is to impose that $\frac{\partial}{\partial t}p_{s,k}[s_{n}]=0$ whenever $x_{s,k}=0$, since there are no nodes of degree $k$ in state $s$.

III Lumping Equations

In this section we present aggregated equations for DBMF and PA. We will first explain all assumptions that are necessary for the lumping, when we partition variables into bins obtained by aggregating degrees according to the degree distribution. We will discuss how to obtain a good partition of the degree set $K$ in Section IV, together with experimental results on several case studies. An important issue of the proposed method is that our aim is to approximate global averages such as the total fraction of nodes that are in a certain state in a given time interval, and not more detailed values such as the probability of (a set of) certain global states of the underlying Markov chain $X$.

III.1 Lumping of Degree-Based Mean-Field Equations

Our next goal is to derive a set of differential equations based on an aggregation of the DBMF equations (see Eq. (II.2)) in order to approximate the total population fraction for each state $s$. The starting point is a partition $B=\{b_{1},\ldots,b_{h}\}\subset 2^{K}$ of the degrees $K=\{1,\ldots,k_{\mathrm{max}}\}$, i.e., such that the disjoint union $\dot{\bigcup}_{b\in B}b$ yields $K$ and $\emptyset\not\in B$. A set of degrees $b\in B$ will be called a bin. Let $x_{s,b}$ denote the corresponding fraction of nodes in state $s$ that have a degree in bin $b$ (also called $(s,b)$ nodes in the following). Knowing the fractions $x_{s,b}$ for all $b\in B,$ we can simply compute $x_{s}=\sum_{b\in B}x_{s,b}P(b)$ to get the total fraction of nodes in state $s$. Here, $P(b)=\sum_{k\in b}P(k)$ is the degree probability of bin $b$. Note that if a network is finite and has $N$ nodes, we can always write $x_{s,b}=\frac{X_{s,b}}{N_{b}},$ where $X_{s,b}$ is the number of $(s,b)$ nodes and $N_{b}$ is the total number of nodes that have degree $k\in b$, i.e., $P(b)=N_{b}/N$. We remark that $x_{s,b}=\sum_{k\in b}x_{s,k}P(k|b)$.

In the next section, we will describe a heuristic to determine a suitable partition $B$. For now, we assume that $B$ is given and that $P(b)\neq 0$ for all $b\in B$. In this case, the probability $P(k|b)=P(k)/P(b)$ that a random node that belongs to bin $b$ is of degree $k$ is well defined. Multiplying both sides of Eq. (6) with $P(k|b)$ and summing over all degrees $k\in b$, we get the following equation for $s\in S$ and $b\in B$.

[TABLE]

Next we replace sums of the form $\sum_{k\in b}x_{s,k}P(k|b)$ by the lumped fraction $x_{s,b}$. However, exact lumping is not possible for terms of the form $\sum_{k\in b}x_{s,k}kP(k|b)$. In order to express the last two sums in terms of the lumped fractions, we make the following assumption:

[TABLE]

In other words, we assume that in bin $b$ the fraction of nodes in state $s$ is the same for each degree $k\in b$. Now, (III.1) becomes

[TABLE]

where $\left\langle k\right\rangle_{b}=\sum_{k\in b}kP(k|b)$ is the average degree in bin $b$. Since the probability $p[s_{n}]$ depends on the fractions $x_{s,k}$, we must also express it in terms of lumped fractions. Imposing Assumption (10) to (3) yields

[TABLE]

Equations (III.1) and (12) define the aggregated equations of the DBMF model.

Since Assumption (10) is usually not true, we obtain approximations of the true fractions $x_{s,b}$ and we use them as approximations of $x_{s,k}$.

III.2 Lumping of Pair Approximation Equations

For the PA equations (6) and (II.3), we also fix a partition $B$ of degrees into bins and multiply Eq. (6) with $P(k|b).$ Summing over all degrees $k\in b$, we get the following lumped equation for the fraction of nodes of bin $b$ that are in state $s.$

[TABLE]

where ${\mu}_{{r_{j_{1}},k}}(r_{j_{2}})$ is again the corresponding first moment of the multinomial distribution (cf.​ Section II.3). Similarly to the DBMF lumping approach, it is necessary to make simplifying assumptions in order to write the last two sums of (III.2) in terms of lumped fractions. More specifically, we again impose Assumption (10) and in addition we assume that the probability that an $(s,k)$ node has an $s_{n}$ neighbor is equal for all $k\in b$, i.e., for all $k\in b$, $p_{s,b}[s_{n}]$ is defined as the probability that a randomly selected neighbor of an $(s,k)$ node is in state $s_{n}\in S$. Hence, if ${\mu}_{s,b}(s_{n})$ is the average number of $s_{n}$ neighbors of an $s$ node in bin $b$, and if this number is binomially distributed, we get

[TABLE]

Consequently, (III.2) gives for all $s\in S$ and $b\in B$,

[TABLE]

Finally, using (14) we get

[TABLE]

To derive a differential equation for the time-dependent variable $p_{s,b}[s_{n}]$ we multiply Equation (II.3) with $P(k|b)$ and sum over all $k\in b$. Assuming $x_{s,k}>0$, this yields

[TABLE]

Using (14) as well as the above assumption $p_{s,k}[s_{n}]=p_{s,b}[s_{n}]$ together with Assumption (10) we get:

[TABLE]

Let $\left\langle k-1\right\rangle_{b}=\sum_{k\in b}(k-1)P(k|b)$. If $r_{j_{2}}\neq s_{n},$ then (see also Appendix D)

[TABLE]

and if $r_{j_{2}}=s_{n}$, then

[TABLE]

We remark that in Eq. (III.2), the calculation of the average rates $\beta^{ss^{\prime}\rightarrow ss_{n}}$ and $\beta^{ss_{n}\rightarrow ss^{\prime}}$ requires the detailed fractions $x_{s,k}$ (cf.​ Eq. (II.3)). Hence, we define, similar to Eqs. (14) and (II.3), the second moments

[TABLE]

and the average rate

[TABLE]

Now, we can replace all unknown terms in (III.2) by using (19) and (20), as well as $\beta^{ss^{\prime}\rightarrow ss^{\prime\prime}}=\beta^{ss^{\prime}\rightarrow ss^{\prime\prime}}_{B}$, where the latter equality holds under Assumption (10), and then the lumped ODE system for the PA model is given by Eqs. (III.2) and (III.2).

IV Binning Method and Results

In this section, we discuss strategies to choose the partition $B$, which strongly influences the quality of the approximation of the lumped equations. We focus on an accurate approximation of the fraction of nodes in state $s$.

Note that if we solve the original DBMF or PA equations we have $x_{s}(t)=\sum_{k\in K}x_{s,k}(t)P(k)$, while from the lumped ODEs we get $\bar{x}_{s}(t)=\sum_{b\in B}\sum_{k\in b}x_{s,b}(t)P(k)$. Hence, the error that we obtain at a certain time $t$ by solving the lumped ODEs instead of the full equations is

[TABLE]

This leads to the following observations.

• •

If $P(k)$ is large, then the deviation $\epsilon_{k,b}(t)$ of $x_{s,b}(t)$ from $x_{s,k}(t)$ will have a great impact on the total error. On the other hand, if $P(k)$ is small, then one could handle even large deviations from the true solution.

• •

For a fixed bin, the smaller the relative difference between the degrees (a more homogeneous bin), the smaller is the error $\epsilon_{k,b}(t)$. For instance, we can expect that the fractions $x_{s,500},x_{s,501},x_{s,502}$ behave much more similarly than the fractions $x_{s,1},x_{s,2},x_{s,3}$ over time. Intuitively, this comes from the fact that for nodes with only one or two neighbors the chance of triggering a contact rule is less similar than the chance of nodes with 500 or 501 neighbors. In other words, the relative difference of the time derivative for degrees $k$ and $k+1$ is smaller for larger degrees than for small ones.

For a given partition $B$, we define the total error $\epsilon_{\text{tot},B}=\max_{t\leq T_{max}}\epsilon_{B}(t),$ where $T_{max}$ is the maximum time horizon.

IV.1 Hierarchical Clustering

We propose a hierarchical clustering approach that takes into account the two observations above and tries to “uniformly” distribute the total error in each bin. Hence, we cluster together consecutive degrees based on their maximum relative distance and on the sum of their probabilities. The latter should be as close as possible to the threshold $1/n,$ where $n$ is the number of bins in the partition. This ensures that no bin contributes disproportionately to the error due to a large probability of the degrees in the bin.

In a bottom-up fashion, we initially consider each degree to form an individual bin and iteratively join those two bins having the smallest distance (as defined below) until a single bin that contains all degrees is obtained. Hence, as a result we get a partition $B_{n}$ for all possible numbers of bins $n$. Note that this computation is very cheap compared to the solution of the ODEs for a fixed partition. Typically, in hierarchical clustering, measures for the distance between two clusters (bins) are based on distance measures between individual points (degrees) and on the choice of a linkage type. In the following, we propose a distance measure that can be directly applied to the bins and takes into account the degree distribution and the bin degrees’ homogeneity. To avoid bins of zero probability (all degrees in the bin have probability zero), we assume that the distance of two bins $b_{i}$ and $b_{j}$, $i\neq j$ is zero if the probability of either of the two bins is zero. Hence, those bins will be merged with other bins. Otherwise, we define the distance $\text{dist}(b_{i},b_{j})$ between two bins $b_{i}$ and $b_{j}$ as a linear combination of their probability and a homogeneity measure of the union $b_{i}\cup b_{j}$, i.e.,

[TABLE]

where $\text{dist}_{hmg}(b_{i},b_{j})$ is a homogeneity measure. Here, $\bar{k}_{b_{i}}$ is the arithmetic mean of bin $b_{i}$, that is $\bar{k}_{b_{i}}=\frac{1}{|b_{i}|}\sum_{k\in b_{i}}k$. Hence the formula measures the distance of the means of two bins after normalizing it with the average of the two means. This compensates for the effect of large differences caused by bins containing large degrees. Finally, $\alpha\in\mathbb{R}_{\geq 0}$ is used to scale $\text{dist}_{hmg}$, since typically $\text{dist}_{hmg}$ yields larger values than $P(b_{i})+P(b_{j})$. We empirically set $\alpha=0.2$; however, the resulting binning turns out to be quite robust to small variations of $\alpha$.

After computing partitions $B_{n}$ for every possible bin number $n$, we choose an adequate $n$ giving a good tradeoff between accuracy and computational time for the solution of the lumped equations. Since we cannot measure the error $\epsilon(n):=\epsilon_{tot,B_{n}}$ between the lumped and the original equations directly, we resort to a simple heuristic to identify a suitable number of bins. This heuristic is based on the function

[TABLE]

where $\mathbf{x}(t,n)$, $\mathbf{x}(t,n-j^{*})$ are the vectors containing the fraction of nodes in each state at time $t$, for the partitions $B_{n},B_{n-j^{*}}$, respectively. The step size ${j}^{*}\geq 1$ allows us to regulate how many steps (bins) are between two solutions we compare. In our implementation, we set ${j}^{*}=5$, to obtain a smoother error function and to more easily detect the convergence of $F(n)$. The basic idea of the heuristic is to increase $n$ and monitor the function $F(n)$, stopping when we observe no noticable variation of its value. Although the partitions that we compute are refinements of each other for increasing $n$, we cannot guarantee that $\epsilon(n)$ and $F(n)$ are always monotonically decreasing. However, $\epsilon(n)$ eventually approaches zero and for all examples that we considered, small changes of $F(n)$ for a number of steps indicated that with refinements of the corresponding partition $B_{n}$, only very small improvements were possible in terms of $\epsilon(n)$. In Fig. 3 (b) and (d) we see that $F(n)$ is roughly proportional to $\epsilon(n)$ and hence the difference $F(n)-F(n-j)$ can be used to identify the number $n$ beyond which $\epsilon(n)$ decreases only very slowly. This proportionality between $F(n)$ and $\epsilon(n)$ was consistent for all models we have studied (see Fig. 3 (b) and (d)).

The above observation gives rise to the heuristic described in Algorithm 1 to decide the number of bins for the lumped equations based on a gradient descent approach. In the algorithm, $j\in\{1,\ldots,n-1\}$ is a fixed (limited) lookahead used to approximate the gradient of $F(n)$, and it is not necessarily equal to the step gap ${j}^{*}$ used in the computation of $F(n)$, even if in our implementation we used the same value, $j=5$, since smaller values (e.g. $j=1$) might not be as effective in indicating convergence. Apart from the model, the network degree distribution and the lookahead $j$, Algorithm 1 takes as input a predefined threshold $\delta$ and a parameter $\gamma_{n}$ which is used for updating the number of bins for the next iteration. Note that increasing $n$ by a fixed number of bins in each step may lead to unnecessarily long running times. When $n$ is small, an increase of the number of bins often gives a large improvement, and hence $-\nabla F(n)$ is large as well. In that case, a large increase of $n$ is adequate. In general, $\gamma_{n}$ can be either fixed for every iteration, i.e., $\gamma_{n}=\gamma$, or it can adaptively be chosen by a backtracking line search method Wright and Nocedal (1999). In our implementation, we fixed $\delta=10^{-3}$ and heuristically chose a fixed scaling constant $\gamma=2\cdot 10^{3}.$

IV.2 Results

We test our lumping approach for two well-studied models, namely the SIR model and a rumor spreading scenario, as well as for a third more complex system. For the simple SIR model, we tested our method for different degree distributions and investigated how the running time of the method increases when we vary the maximum degree $k_{\mathrm{max}}$ of the network. For the rumor spreading model, we show that even with a small number of bins the method is able to capture contact processes with more than one contact rules. Finally, using the complex disease model, we verify that the method performs well even in cases of networks with many possible node states. For all models, we compare the dynamics and the running times of the original equations versus the lumped equations for a number of bins needed to get an error of at most $3\cdot 10^{-3}$. Note that for each model we chose, whenever possible, parameter values s.t. the original DBMF equations differ significantly from the PA equations. In addition, we report the number of bins that was returned considering the stopping criterion of Algorithm 1. In well, we provide a user-friendly python tool, called LUMPY LUMPY Tool (2017), that takes a simple description of a contact process as input and a degree distribution. It automatically generates all kinds of equations, both the full ones and the lumped ones, using the binning heuristic discussed above.

SIR Model

The first model we examined is the SIR model, which consists of the following rules: a susceptible node $S$ that interacts with an infected node gets infected at rate $\theta_{1}$; an infected node can recover at rate $\theta_{2}$; and a recovered node can go back to the susceptible state at rate $\theta_{3}$:

[TABLE]

Here for all the experiments we set $\theta=(\theta_{1},\theta_{2},\theta_{3})=(6,4,1)$ and we assume the underlying network’s structure follows a power law distribution with $k_{\mathrm{max}}=10^{3}$ and $P(k)\propto k^{-2.4}$. In Fig. 3 we show the comparison between the dynamics of the aggregated DBMF equations using 20 bins and the original DBMF equations and in Fig. 3 we compare the aggregated PA equations for 21 bins and the original PA equations. We see that it in both cases is almost impossible to distinguish between the global measures derived from the aggregated ODEs and those obtained from the original system of equations. In Fig. 3 and 3 we plot the maximum error $\epsilon(n)$ between the full and the lumped DBMF and PA equations, respectively, as a function of the number of bins $n$ as well as $F(n)$, i.e., the maximum difference between the aggregated solutions for $n$ and $n-j$ ($j=5$). Note that both errors quickly converge to zero as the number of bins increases.

In Fig. 4 and 4 we show the same plots for a network with uniform degree distribution $P(k)=1/k_{\mathrm{max}},$ where $k_{\mathrm{max}}=500$, and the solution of the PA equations (aggregated vs. original ones). For the chosen network size the dynamics of DBMF equations completely agree with those of PA. Even though the uniform degree distribution is unrealistic for a typical contact network, we see that the lumped equations give very accurate results for a small number of bins. For an error $\epsilon(n)$ of at most $3\cdot 10^{-3}$ we only need $n=11$ bins and $132$ ODEs, in comparison to $500$ degrees and $6000$ ODEs of the full PA equations for this network. Note that the approximation works that accurately here, although due to the uniform assumption only the second part of the distance function in Eq. (23) (homogeneity of degrees) is considered in clustering.

In Fig. 5 we plot the number of bins needed for the solution of the DBMF and PA equations for an approximation error of at most $3\cdot 10^{-3}$ in a power law network with varying $k_{\mathrm{max}}.$ The network size grows exponentially, while the number of bins needed grows only linearly. It is worth mentioning that in order to be sure that the dynamics of the network will change as $k_{\mathrm{max}}$ increases, we have changed the parameter $\alpha$ of the power law in such a way that the coefficient of variation, i.e., $\sigma/\!\left\langle k\right\rangle$, of the distribution remains constant for different $k_{\mathrm{max}}$ (see also Appendix, E).

In Fig. 5 we plot the gain ratio $T_{full}/T_{lumped}$ of the computation times for solving the aggregated equations instead of the full system. Generally the speedups are impressively high, and we see a significant increase of the ratio as the network size increases. In particular, for $k_{\mathrm{max}}=10^{6}$ we need $5.2$ days for solving the full equations, versus $12$ sec.​ for the lumped equations. The precise computation times of the needed lumped equations for varying $k_{\mathrm{max}}$ are listed in Table 1, where the additional time needed for the hierarchical clustering is in the order of seconds up to network size $k_{\mathrm{max}}=10^{4}$ and in the order of minutes up to $k_{\mathrm{max}}=10^{6}.$

In Table 2 we include the number of bins returned by Algorithm 1, together with the corresponding error, for both DBMF and PA equations for varying $k_{\mathrm{max}}$. The number of bins returned by Algorithm 1 for $\delta=10^{-3}$ yields very accurate approximations (error $\ll$ $3\cdot 10^{-3}$), as the chosen stopping criterion of the heuristic is rather conservative. In any case, the speedups compared to a full solution of the equations are very similar to those in Fig. 5.

Rumor Spreading Model

As a second model, we consider a rumor spreading model Barrat et al. (2008) that includes three different types of individuals: ignorants (I), spreaders (S), and stiflers (R). An ignorant that is in contact with a spreader becomes a spreader at rate $\lambda.$ A spreader becomes a stifler when the spreader is in contact either with a spreader or a stifler; in both cases the corresponding transition rate is $\alpha$:

[TABLE]

Here, we set $\lambda=6.0$ and $\alpha=0.7$ and, for the network’s degree distribution, we consider a truncated power law $P(k)\propto k^{-2.2}$ with $k_{\mathrm{max}}=500.$ In Fig. 6 we show the dynamics of the full DBMF ODEs compared to the DBMF lumped equations and in Fig. 6 the corresponding comparison for the PA equations. We note that DBMF here gives a poor approximation, which is qualitatively different from the PA solution. Even though there is a great difference between the solutions of the two systems of equations, the lumped equations approximate very accurately the original equations for both cases. For this model, we needed only 13 bins and 39 ODEs for the DBMF case and 27 bins and 324 ODEs for the PA in order to get a maximum error of $3\cdot 10^{-3}$. The heuristic in Alg. 1 with $\delta=10^{-3}$ returned 43 bins with an error of $10^{-4}$ for the first and 60 bins with an error of $2\cdot 10^{-4}$ for the second case. The computation times for solving the lumped equations of DBMF and PA for every possible number of bins are shown in Fig. 6 (b) and (d). Again, we see huge computational gains when solving the lumped equations instead, especially in the case of PA.

SIIIR Model

Recent infection models describe the infection process in more detail by considering multiple infection phases Miller and Volz (2013). In HIV, for instance, the early phase is well known to be significantly more infectious than a longer-lasting chronic phase Pilcher et al. (2004); Eames and Keeling (2002). In order to capture such behavior, here we incorporate the phases as different states of the individuals, i.e., $I,II,III,$ with different transmission rates:

[TABLE]

We fix $\theta=(\theta_{1},\ldots,\theta_{7})=(5,1.5,1,2,2,2,2)$ and we test our lumped equations with the original DBMF and PA solutions taken for the scale-free network $(P(k)\propto k^{-2.5},k_{\mathrm{max}}=100)$ of sexual partners over a time period of one year Liljeros et al. (2001). Despite the large number of possible node states, the dynamics of the total fractions of all states is captured very accurately (error of less than $3\cdot 10^{-3}$) by the lumped equations with 20 bins in both DBMF (Fig. 7) and PA (Fig. 7). For this system, Algorithm 1 returned 39 bins for DBMF and 46 bins for PA with an error of $4.7\cdot 10^{-5}$ and $5\cdot 10^{-5}$, respectively, when $\delta=10^{-3}$ was chosen. Last, it is worth mentioning that we also tested the robustness of our method regarding the choice of the initial fractions. In Gleeson (2011) and Gleeson (2013) the author uses for each state $s$ a constant initial fraction for all $x_{s,k}$ for every $k$, which is equal to the initial total fraction of $x_{s}.$ Here, we considered in addition completely random initial fractions for each $x_{s,k}$ for all states $s$, and this initialization returned an even smaller error function in the beginning (see also Appendix, F), mainly because we avoid the extreme values that might be present in our initialization, such as the total fraction of susceptible nodes, which is here close to one. Hence, we conclude that the method is robust against different initial conditions.

V Conclusion

Studying complex dynamic processes happening on networks is a daunting computational task when relying on simulation. To ameliorate this problem, several forms of approximations have been considered, all trying to take the topology of the network into account to some extent: degree-based mean-field Pastor-Satorras and Vespignani (2001), AME, and pair approximation Gleeson (2013). AME and PA equations, in particular, seem to capture with precision the dynamics of a large variety of network topologies, but in Gleeson (2013) they are typically restricted to the binary case, i.e., to processes with only two local states. Here, we consider the AME and the subsequent derivation of the PA in the general case of contact processes of $n$ states. The total number of equations for both DBMF and PA is directly dependent on the maximum degree, $|S|k_{\mathrm{max}}$ and $|S|^{2}k_{\mathrm{max}}$, respectively, and hence the total computation time even for these approximate methods can be demanding for large networks with many states.

In this paper, we argue that although the network structure is proven to be essential for its dynamics, we do not need a full resolution to get precise results for the time evolution of the total fractions. For instance, it is expected that the degrees belonging to the tail of a power law should not significantly change the precision of the equations if clustered together. Motivated by this insight, we investigated whether it is possible to capture the dynamics of DBMF and PA equations by solving a smaller set of equations. Taking advantage of the degree distribution characteristics, and additionally exploiting a notion of homogeneity of the differences between separate degrees, we provide a mechanism according to which we can effectively reduce the resolution of network degrees, clustering them into a small number of bins and constructing a small set of aggregated differential equations.

The results in terms of computational gain are particularly impressive: we manage to approximate collective network statistics, namely the total fractions of nodes in each state, produced by DBMF and PA equations by solving a significantly smaller number of differential equations and gaining a computational speedup of several orders of magnitude. We are able to reproduce these results for a number of different network degree distributions and a number of different models.

Furthermore, the procedure to construct the reduced set of equations is fully automatic, and it is implemented as an open-source Python tool LUMPY Tool (2017). The tool takes a compact description of the model as input and returns Python files containing the lumped and the full DBMF and PA equations as output. This makes validation of the results and experimentation with other models a very easy task.

The natural direction in which to extend this work is to provide a similar approach to aggregate the AME equations. A direct strategy would be to use the same binning as for PA equations and lump AME only over different degrees $k$. For each of the formed bins, we would then need to solve $|S||\mathbb{M}_{k^{b}_{\mathrm{max}}}|$ equations. This is guaranteed to give at least as good results as the lumped PA equations. The quality, however, may depend on whether the solution of the original PA equations is close to the solution of the AME. However, getting a simultaneous aggregation over degrees $k\in K$ and neighbors $\mathbb{m}\in\mathbb{M}_{k}$ which is close to the solution of full AME equations could give a tremendous computational boost and make AME tractable regardless of the network size. Finally, a modification of the presented lumping method could successfully be applied to adaptive networks, where the degree distribution changes smoothly over time. In this case, we expect that “on-the-fly” and local rearrangements of the clusters could be sufficient to capture well the topology changes of such networks.

Acknowledgements.

We would like to thank Mayank Goyal for an early implementation of the binning approach on the binary state case dynamics as part of his master’s thesis.

Appendix A State distribution of the neighborhood in uncorrelated networks

In this section, we prove Eq. (3) for uncorrelated networks, i.e., $p_{k}[s^{\prime}]=p[s^{\prime}]=1/{\left\langle k\right\rangle}\sum_{k^{\prime}\in K}{P(k^{\prime})k^{\prime}x_{s^{\prime},k^{\prime}}}$.

We consider first the probability $q_{k}(k^{\prime})$ that, given a node of degree $k$, a randomly chosen direct neighbor of it has degree $k^{\prime}$. Since in an uncorrelated network the above probability is independent of $k$, it is the same as the probability $q(k^{\prime})$ of choosing a random edge that connects to a node of degree $k^{\prime}$. Hence:

[TABLE]

where $N_{k^{\prime}}$ is the number of nodes of degree $k^{\prime}$ and $|E|$ is the total number of edges. Please note here the subtle but important difference between $q(k^{\prime})$, which is the probability of picking a $k^{\prime}$ node after randomly choosing a network edge, and $P(k^{\prime})$, which is the probability of a random node being of degree $k^{\prime}$. Thus the probability that a random neighbor of a node of degree $k$ in an uncorrelated network is in state $s^{\prime}$ is independent of $k$, and given by:

[TABLE]

Appendix B Number of AME equations

Without considering the redundant equations, the number of equations of the general AME case with $|S|$ states is $\binom{k_{\mathrm{max}}+|S|}{|S|-1}(k_{\mathrm{max}}+1).$

Proof:

[TABLE]

The first equality is the rising sum of binomial coefficients and the rest is simple calculus.

Appendix C From AME to PA

Proof: We consider the general case of the AME equation in (4) for $|S|$ states and $N_{I}$ + $N_{C}$ rules. Assuming that the number of links of an $(s,k)$ node to all other states is multinomially distributed, we get two equations: The first equation describes the fraction of $(s,k)$ nodes and the second equation the probability $p_{s,k}[s_{n}]$ that a random neighbor of $(s,k)$ is in state $s_{n}$.

$\circ$ In order to derive an expression for the fraction of $(s,k)$ nodes, we sum up both parts of (4) over all possible vectors $\mathbb{m}\in\mathbb{M}_{k},$ where $\mathbb{M}_{k}$ is the support of the multinomial distribution, and we apply the multinomial assumption to all terms. This yields

[TABLE]

which after simplification becomes

[TABLE]

Here, ${\mu}_{s,k}(s_{i})$ denotes the average of the $i$-th entry over all $\mathbb{m}$ vectors of $(s,k)$ nodes. Applying the multinomial assumption, it holds that ${\mu}_{s,k}(s_{i})=\sum_{\mathbb{m}\in\mathbb{M}_{k}}f(\mathbb{m},p_{s,k})\mathbb{m}[s_{i}]=kp_{s,k}[s_{i}],$ where $f(\mathbb{m},p_{s,k})$ is the probability of picking a neighbor vector $\mathbb{m}$ from a multinomial distribution with a vector of success probabilities equal to $p_{s,k}.$

Note that the last two terms in (C) cancel out because

[TABLE]

where the last equality follows from the fact that when $\mathbb{m}[s^{\prime\prime}]=k$ then $\mathbb{m}[s^{\prime}]=0$.

$\circ$ To derive an equation for $p_{s,k}[s_{n}]$, we first consider the change over time of the number of edges between a node in state $s$ and $s_{n}$. From the multinomial assumption, we have that

[TABLE]

where $f(\mathbb{m},p_{s,k})$ is the probability of picking a neighbor vector $\mathbb{m}$ from a multinomial distribution with a vector of success probabilities equal to $p_{s,k}.$ Using similar arguments, we can now derive an expression for (29) by summation of (4) after multiplication of $\mathbb{m}[s_{n}]$.

[TABLE]

For the last two sums, we distinguish the cases (i) $s_{n}\neq s^{\prime}\neq s^{\prime\prime}$ and (ii) $s_{n}=s^{\prime}\vee s_{n}=s^{\prime\prime}.$

In case (i) we get

[TABLE]

which implies that the difference (1)-(2) is zero for this case. Note that the entry of $s_{n}$ in vector $\mathbb{m}$ remains the same during substitution, since $s_{n}\neq s^{\prime}\neq s^{\prime\prime}$.

In case (ii) we distinguish again between the two subcases: (a) $s_{n}=s^{\prime}$, and (b) $s_{n}=s^{\prime\prime}.$

Subcase (a): The difference (1)-(2) becomes

[TABLE]

Subcase (b): Similar to subcase (a), the difference (1)-(2) becomes

[TABLE]

Thus, overall, dividing both parts of (C) by $x_{s,k}k$, rearranging the term and doing all the possible algebraic simplifications, we get

[TABLE]

Next we replace and get:

[TABLE]

Finally, after imposing the multinomial assumption in the computation of the rates $\beta^{ss^{\prime}\rightarrow ss^{\prime\prime}}$, we get

[TABLE]

For $s,s_{n},s^{\prime}_{n}\in S$ and $k\in K$, we have the following:

• •

The expectation of the multinomial distribution is

[TABLE]

• •

For the second mixed moment we distinguish the cases

1. (i)

$s_{n}\neq s^{\prime}_{n}$

[TABLE] 2. (ii)

$s_{n}=s^{\prime}_{n}$

[TABLE]

In the above equations $\bar{\mu}_{s,k}(s_{n},s_{n}^{\prime})$ is the second central mixed moment, i.e., the covariance of the multinomial distribution for the number of edges between an $(s,k)$ node and $s_{n}$ nodes and between an $(s,k)$ node and $s_{n}^{\prime}$ nodes, while $\bar{\mu}_{s,k}(s_{n},s_{n})$ is the second central moment, i.e., the variance for the number of edges between an $(s,k)$ node and $s_{n}$ nodes.

Appendix D Simplification of sums for lumping PA

To simplify the sum $\sum_{k\in b}\frac{1}{k}{\mu}_{s,k}(r_{j_{2}},s_{n})P(k|b)$ of Eq. (III.2), we consider the two cases of (37) and (38) for the second mixed moment and we get:

1. (i)

$s_{n}\neq s^{\prime}_{n}$

[TABLE] 2. (ii)

$s_{n}=s^{\prime}_{n}$

[TABLE]

Similarly, for simplifying ${\mu}_{s,b}(s_{n},s_{n}^{\prime})$ of Eq. (21) that is used for the computation of $\beta^{ss^{\prime}\rightarrow ss^{\prime\prime}}_{B}$, we get:

1. (i)

$s_{n}\neq s^{\prime}_{n}$

[TABLE] 2. (ii)

$s_{n}=s^{\prime}_{n}$

[TABLE]

Appendix E DBMF dynamics for SIR model

In Fig. 9 we show the dynamics of the SIR models in several networks with $k_{\mathrm{max}}\in\{10^{2},10^{3},10^{5},10^{6}\}$. The degree distribution of each of the networks was chosen such that the coefficient of variation $\sigma/\left\langle k\right\rangle$ of the distribution remains constant. For a truncated power law probability distribution $P(k)=ck^{-\alpha},k\in[k_{\mathrm{min}},k_{\mathrm{max}}]$ where $c$ is the normalizing constant $c=\frac{\alpha-1}{k_{\mathrm{min}}^{-\alpha+1}},$ it holds that \left\langle k\right\rangle=\alpha\Big{[}\frac{k^{2-\alpha}}{2-\alpha}\Big{]}^{k_{\mathrm{max}}}_{k_{\mathrm{min}}} and \left\langle k^{2}\right\rangle=\alpha\Big{[}\frac{k^{3-\alpha}}{3-\alpha}\Big{]}^{k_{\mathrm{max}}}_{k_{\mathrm{min}}}. We get $\sigma=\sqrt{\left\langle k^{2}\right\rangle-\left\langle k\right\rangle{{}^{2}}},$ and fixing $\sigma/\left\langle k\right\rangle=7.5$ and $k_{\mathrm{min}}=1$, we compute $\alpha$ for each different $k_{\mathrm{max}}.$

Appendix F Initial conditions for SIIIR model

In Fig. 9 we show the PA lumping error of the SIIIR model in a power law network with $k_{\mathrm{max}}=100,P(k)\propto k^{-2.5}$ over time for different choices of the initial fractions $x_{s,k},$ where $s\in\{S,I,II,III,R\}$ and $k\in K$. We plot the error of the initial conditions we used for the results of Fig. 7 (our initial conditions), the error of an initialization that gives the same values to all $x_{s,k}$ for all $s$ and $k$ (common uniform), and finally we plot the average error of 30 different random initial conditions for the fractions $x_{s,k}$ together with the error’s sample standard deviation $(\pm\sigma).$