Reducing Spreading Processes on Networks to Markov Population Models

TL;DR

This paper introduces a novel lumping scheme that reduces complex network-based epidemic models to Markov Population Models, enabling more efficient analysis and approximation of spreading processes.

Contribution

The authors propose a new node-partitioning lumping method that transforms complex network epidemic models into Markov Population Models, facilitating the use of existing approximation techniques.

Findings

Lumping reduces the state space size significantly.

Different counting abstractions affect approximation accuracy.

Numerical examples demonstrate the method's effectiveness.

Abstract

Stochastic processes on complex networks, where each node is in one of several compartments, and neighboring nodes interact with each other, can be used to describe a variety of real-world spreading phenomena. However, computational analysis of such processes is hindered by the enormous size of their underlying state space. In this work, we demonstrate that lumping can be used to reduce any epidemic model to a Markov Population Model (MPM). Therefore, we propose a novel lumping scheme based on a partitioning of the nodes. By imposing different types of counting abstractions, we obtain coarse-grained Markov models with a natural MPM representation that approximate the original systems. This makes it possible to transfer the rich pool of approximation techniques developed for MPMs to the computational analysis of complex networks' dynamics. We present numerical examples to investigate the relationship between the accuracy of the MPMs, the size of the lumped state space, and the type of counting abstraction.

Equations

$$\mathcal{X}=\{x\mid x:\mathcal{N}\to \mathcal{S}\}$$

$$\mathcal{M}=\Big{\{}\mathbf{m}\in\mathbb{Z}_{\geq 0}^{|\mathcal{S}|}\mathrel{\bigg{|}}\sum_{s\in\mathcal{S}}\mathbf{m}[\text{s}]\leq k_{\text{max}}\Big{\}}\;,$$

$$\mathtt{S}\stackrel{f}{\to }\mathtt{I}\text{with}f(\mathbf{m})=\lambda \cdot \mathbf{m}[\mathtt{I}],$$

$$\mathtt{I}\stackrel{f}{\to }\mathtt{S}\text{with}f(\mathbf{m})=\mu ,$$

$$\mathcal{L}:\mathcal{X}\to \mathcal{Y}$$

$$q(y,y^{\prime })=\frac{1}{|{\mathcal{L}}^{-1}(y)|}\sum _{x\in {\mathcal{L}}^{-1}(y)}\sum _{x^{\prime }\in {\mathcal{L}}^{-1}(y^{\prime })}q(x,x^{\prime }).$$

$$\mathcal{Y}=\{y\mid y:\mathcal{S}\times \mathcal{P}\to {\mathbb{Z}}_{\ge 0}\}$$

$$\mathcal{L}(x)=y$$

$$y(s,P)=|\{n\in \mathcal{N}\mid X(n)=s,n\in P\}|.$$

$$\mathcal{N}:=\mathcal{N}\cup \{n_{\star }\}\mathcal{S}:=\mathcal{S}\cup \{\star \}L(n_{\star })=\star \mathcal{P}:=\mathcal{P}\cup \{P_{\star }\}$$

$$\mathcal{E}:=\mathcal{E}\cup \{(n,n_{\star })\mid n\in \mathcal{N},n\ne n_{\star }\}.$$

$$\mathcal{Y}=\{y\mid y:\mathcal{S}\times \mathcal{P}\times \mathcal{S}\times \mathcal{P}\to {\mathbb{Z}}_{\ge 0}\}$$

$$\displaystyle\mathcal{L}(x)=y\phantom{\Big{(}}$$

$$\displaystyle y(s,P,s^{\prime},P^{\prime})=|\big{\{}(n,n^{\prime})\in\mathcal{E}\mid x(n)=s,n\in P,x(n^{\prime})=s^{\prime},n^{\prime}\in P^{\prime}\big{\}}|$$

$$d_{\text{k}}(n,n^{\prime })=\frac{|k_n-k_{n^{\prime }}|}{\text{max}(k_n,k_{n^{\prime }})}.$$

$$\displaystyle q(\mathbf{y},\mathbf{y}^{\prime})=\begin{cases*}\alpha(\mathbf{y})&if $\exists(\alpha,\mathbf{b})\in\mathcal{R},\mathbf{y}^{\prime}=\mathbf{y}+\mathbf{b}$\\ 0&otherwise\end{cases*}\;.$$

$$\mathcal{Z}=\{(s,P)\mid s\in \mathcal{S},P\in \mathcal{P}\}.$$

$${\alpha }_{r,P}:$$

$${\alpha }_{r,P}(\mathbf{y})=$$

$$\displaystyle\mathbf{b}_{r,P}[z]=\begin{cases*}1&if $z.s=s_{2},P=z.P$\\ -1&if $z.s=s_{1},P=z.P$\\ 0&otherwise\end{cases*}\;.$$

$$\displaystyle\mathcal{Z}=\big{\{}(s_{source},P_{source},s_{target},P_{target})\mid$$

$$\displaystyle(s_{source},P_{source})\leq(s_{target},P_{target})\big{\}}\;.$$

$${\mathcal{V}}_P=\bigcup _{n\in P}{\mathcal{V}}_n.$$

$${\alpha }_{r,P,\mathbf{v}}:$$

$${\alpha }_{r,P,\mathbf{v}}(\mathbf{y})=$$

$$\displaystyle\mathbf{b}_{r,P,\mathbf{v}}[z]=\begin{cases*}\mathbf{v}[z.s_{target},z.P_{target}]&if $s_{2}=z.s_{source},P=z.P_{source}$\\ -\mathbf{v}[z.s_{target},z.P_{target}]&if $s_{1}=z.s_{source},P=z.P_{source}$\\ \mathbf{v}[z.s_{source},z.P_{source}]&if $s_{2}=z.s_{target},P=z.P_{target}$\\ -\mathbf{v}[z.s_{source},z.P_{source}]&if $s_{1}=z.s_{target},P=z.P_{target}$\\ 0&otherwise\end{cases*}\;.$$

$$|\mathcal{Y}|=\text{scaleobj}1.2\prod _{P\in \mathcal{P}}\left(\genfrac{}{}{0ex}{}{|P|+|\mathcal{S}|-1}{|\mathcal{S}|-1}\right).$$

$$|\mathcal{Y}|\le \text{scaleobj}1.1\prod _{\begin{array}{c}{P}^{\prime },P^{\prime \prime }\in {\mathcal{P}}^2\\ P^{\prime }\le P^{\prime \prime }\end{array}}\left(\genfrac{}{}{0ex}{}{ϵ(P^{\prime },P^{\prime \prime })+{\mathcal{S}}^2-1}{{\mathcal{S}}^2-1}\right)\cdot \text{scaleobj}1.1\prod _{P\in \mathcal{P}}\left(\genfrac{}{}{0ex}{}{|P|+|\mathcal{S}|-1}{|\mathcal{S}|-1}\right).$$

$$P_{L}\big{(}Y(t)=x\big{)}=\frac{P\big{(}Y(t)=y\big{)}}{|\mathcal{L}^{-1}(y)|}\text{\hskip 14.22636pt}\text{where $y$ is s.t.\;$L(x)=y$.}$$

$$d(P,P_{L})=\max_{t}\sum_{x\in\mathcal{X}}\Big{|}P_{L}\big{(}Y(t)=x)-P(X(t)=x\big{)}\Big{|}\;.$$

$$\displaystyle\alpha_{r,P}(\mathbf{y})=\sum\limits_{n\in P}\phantom{.}\sum\limits_{\mathbf{v}\in\mathcal{V}_{n}}f(\mathbf{m}_{\mathbf{v}})\Pr\big{(}X(n)=s_{1},V(n)=\mathbf{v}\big{)}\;,$$

$$\displaystyle\Pr\Big{(}X(n)=s_{1},V(n)=\mathbf{v}\Big{)}$$

$$=$$

$$\displaystyle\Pr\Big{(}X(n)=s_{1}\Big{)}=\frac{\mathbf{y}[s_{1},P]}{|P|}\hskip 14.22636pt\text{where: }n\in P\;.$$

$${\mathbf{y}}_P[s]=\mathbf{y}[s,P].$$

Full text

In this work, we demonstrate that lumping can be used to reduce any epidemic model to a Markov Population Model (MPM). Therefore, we propose a novel lumping scheme based on a partitioning of the nodes. By imposing different types of counting abstractions, we obtain coarse-grained Markov models with a natural MPM representation that approximate the original systems. This makes it possible to transfer the rich pool of approximation techniques developed for MPMs to the computational analysis of complex networks’ dynamics.

We present numerical examples to investigate the relationship between the accuracy of the MPMs, the size of the lumped state space, and the type of counting abstraction.

1 Introduction

Computational modeling and analysis of dynamic processes on networked systems is a wide-spread and thriving research area. In particular, much effort has been put into the study of spreading phenomena [2, 36, 15, 26]. Arguably, the most common formalism for spreading processes is the so-called Susceptible-Infected-Susceptible (SIS) model with its variations [26, 36, 37].

In the SIS model, each node is either infected (I) or susceptible (S). Infected nodes propagate their infection to neighboring susceptible nodes and become susceptible again after a random waiting time. Naturally, one can extend the number of possible node states (or compartments) of a node. For instance, the SIR model introduces an additional recovered state in which nodes are immune to the infection.

SIS-type models are remarkable because—despite their simplicity—they allow the emergence of complex macroscopic phenomena guided by the topological properties of the network. There exists a wide variety of scenarios which can be described using the SIS-type formalism. For instance, the SIS model has been successfully used to study the spread of many different pathogens like influenza [24], dengue fever [38], and SARS [34]. Likewise, SIS-type models have shown to be extremely useful for analyzing and predicting the spread of opinions [47, 27], rumors [51, 50], and memes [49] in online social networks. Other areas of applications include the modeling of neural activity [14], the spread of computer viruses [10] as well as blackouts in financial institutions [32].

The semantics of SIS-type processes can be described using a continuous-time Markov chain (CTMC) [26, 45] (cf. Chapter 3 for details). Each possible assignment of nodes to the two node states S and I constitutes an individual state in the CTMC (here referred to as network state to avoid confusion111In the following, we will use the term CTMC state and network state interchangeably.). Hence, the CTMC state space grows exponentially with the number of nodes, which renders the numeral solution of the CTMC infeasible for most realistic contact networks.

This work investigates an aggregation scheme that lumps similar network states together and thereby reduces the size of the state space. More precisely, we first partition the nodes of the contact network. After which, we impose a counting abstraction on each partition. We only lump two networks states together when their corresponding counting abstractions coincide on each partition.

As we will see, the counting abstraction induces a natural representation of the lumped CTMC as a Markov Population Model (MPM). In an MPM, the CTMC states are vectors which, for different types of species, count the number of entities of each species. The dynamics can elegantly be represented as species interactions. More importantly, a very rich pool of approximation techniques has been developed on the basis of MPMs, which can now be applied to the lumped model. These include efficient simulation techniques [6, 1], dynamic state space truncation [22, 31], moment-closure approximations [42, 18], linear noise approximation [44, 17], and hybrid approaches [3, 41].

The remainder of this work is organized as follows: Section 2 shortly revises related work, Section 3 formalized SIS-type models and their CTMC semantics. Our lumping scheme is developed in Section 4. In Section 5, we show that the lumped CTMCs have a natural MPM representation. Numerical results are demonstrated in in Section 6 and some conclusions in Section 7 complete the paper and identify open research problems.

2 Related Work

The general idea behind lumping is to reduce the complexity of a system by aggregating (i.e., lumping) individual components of the system together. Lumping is a popular model reduction technique which has been used to reduce the number of equations in a system of ODEs and the number of states in a Markov chain, in particular in the context of biochemical reaction networks [29, 5, 48, 7]. Generally speaking, one can distinguish between exact and approximate lumping [29, 5].

Most work on the lumpability of epidemic models has been done in the context of exact lumping [26, 40, 46]. The general idea is typically to reduce the state space by identifying symmetries in the CTMC which themselves can be found using symmetries (i.e., automorphisms) in the contact network. Those methods, however, are limited in scope because these symmetries are infeasible to find in real-world networks and the state space reduction is not sufficient to make realistic models small enough to be solvable.

This work proposes an approximate lumping scheme. Approximate lumping has been shown to be useful when applied to mean-field approximation approaches of epidemic models like the degree-based mean-field and pair approximation equations [28], as well as the approximate master equation [19, 13]. However, mean-field equations are essentially inflexible as they do not take topological properties into account or make unrealistic independence assumptions between neighboring nodes.

Moreover, [25] proposed using local symmetries in the contact network instead of automorphisms to construct a lumped Markov chain. This scheme seems promising, in particular on larger graphs where automorphisms often do not even exist, however, the limitations for real-world networks due to a limited amount of state space reduction and high computational costs seem to persist.

Conceptually similar to this work is also the unified mean-field framework (UMFF) proposed by Devriendt et al. in [9]. Devriendt et al. also partition the nodes of the contact network but directly derive a mean-field equation from it. In contrast, this work focuses on the analysis of the lumped CTMC and its relation to MPMs. Moreover, we investigate different types of counting abstractions, not only node based ones.

3 Spreading Processes

Let $\mathcal{G}=(\mathcal{N},\mathcal{E})$ be a an undirected graph without self-loops. At each time point $t\in\mathbb{R}_{\geq 0}$ each node occupies one of $m$ different node states, denoted by $\mathcal{S}=\{s_{1},s_{2},\dots,s_{m}\}$ (typically, $\mathcal{S}=\{\mathtt{S},\mathtt{I}\})$. Consequently, the network state is given by a labeling $x:\mathcal{N}\rightarrow\mathcal{S}$. We use

[TABLE]

to denote all possible labelings. $\mathcal{X}$ is also the state space of the underlying CTMC. As each of the $|\mathcal{N}|$ nodes occupies one of $m$ states, we find that $|\mathcal{X}|=|\mathcal{S}|^{|\mathcal{N}|}$.

A set of stochastic rules determines the particular way in which nodes change their corresponding node states. Whether a rule can be applied to a node depends on the state of the node and of its immediate neighborhood.

The neighborhood of a node is modeled as a vector $\mathbf{m}\in\mathbb{Z}_{\geq 0}^{|\mathcal{S}|}$ where $\mathbf{m}[s]$ denotes the number of neighbors in state $s\in\mathcal{S}$ (we assume an implicit enumeration of states). Thus, the degree (number of neighbors, denoted by $k$) of a node is equal to the sum over its associated neighborhood vector, that is, $k=\sum_{s\in\mathcal{S}}\mathbf{m}[s]$. The set of possible neighborhood vectors is denoted as

[TABLE]

where $k_{\text{max}}$ denotes the maximal degree in a given network.

Each rule is a triplet $s_{1}\xrightarrow{f}s_{2}$ ($s_{1},s_{2}\in\mathcal{S},s_{1}\neq s_{2}$), which can be applied to each node in state $s_{1}$. When the rule \sayfires it transforms the node from $s_{1}$ into $s_{2}$. The rate at which a rule \sayfires is specified by the rate function $f:\mathcal{M}\rightarrow\mathbb{R}_{\geq 0}$ and depends on the node’s neighborhood vector. The time delay until the rule is applied to the network state is drawn from an exponential distribution with rate $f(\mathbf{m})$. Hence, higher rates correspond to shorter waiting times. For the sake of simplicity and without loss of generality, we assume that for each pair of states $s_{1}$, $s_{2}$ there exists at most one rule that transforms $s_{1}$ to $s_{2}$.

In the well-known SIS model, infected nodes propagate their infection to susceptible neighbors. Thus, the rate at which a susceptible node becomes infected is proportional to its number of infected neighbors:

[TABLE]

where $\lambda\in\mathbb{R}_{\geq 0}$ is a rule-specific rate constant (called infection rate) and $\mathbf{m}[\mathtt{I}]$ denotes the number of infected neighbors. Furthermore, a recovery rule transforms infected nodes back to being susceptible:

[TABLE]

where $\mu\in\mathbb{R}_{\geq 0}$ is a rule-specific rate constant called recovery rate.

A variation of the SIS model is the SI model where no curing rule exists and all nodes (that are reachable from an infected node) will eventually end up being infected. Intuitively, each rule tries to \sayfire at each position $n\in\mathcal{N}$ where it can be applied. The rule and node that have the shortest waiting time \saywin and the rule is applied there. This process is repeated until some stopping criterion is fulfilled.

3.1 CTMC Semantics

Formally, the semantics of the SIS-type processes can be given in terms of continuous-time Markov Chains (CTMCs). The state space is the set of possible network states $\mathcal{X}$. The CTMC has a transition from state $x$ to $x^{\prime}$ ($x,x^{\prime}\in\mathcal{X}$, $x\neq x^{\prime}$) if there exists a node $n\in\mathcal{N}$ and a rule $s_{1}\xrightarrow{f}s_{2}$ such that the application of the rule to $n$ transforms the network state from $x$ to $x^{\prime}$. The rate of the transition is exactly the rate $f(\mathbf{m})$ of the rule when applied to $n$. We use $q(x,x^{\prime})\in\mathbb{R}_{\geq 0}$ to denote the transition rate between two network states. Fig. 1 illustrates the CTMC corresponding to an SIS process on a small toy network.

Explicitly computing the evolution of the probability of $x\in\mathcal{X}$ over time with an ODE solver, using numerical integration, is only possible for very small contact networks, since the state space grows exponentially with the number of nodes. Alternative approaches include sampling the CTMC, which can be done reasonably efficiently even for comparably large networks [20, 8, 43] but is subject to statistical inaccuracies and is mostly used to estimate global properties.

4 Approximate Lumping

Our lumping scheme is composed of three basic ingredients:

Node Partitioning: The partitioning over the nodes $\mathcal{N}$ that is explicitly provided.

Counting Pattern: The type of features we are counting, i.e., nodes or edges.

Implicit State Space Partitioning: The CTMC state space is implicitly partitioned by counting the nodes or edges on each node partition.

We will start our presentation discussing the partitioning of the state space, then showing how to obtain it from a given node partitioning and counting pattern. To this end, we use $\mathcal{Y}$ to denote the new lumped state space and assume that there is a surjective222If $\mathcal{L}$ is not surjective, we consider only the image of $\mathcal{L}$ to be the lumped state space. lumping function

[TABLE]

that defines which network states will be lumped together. Note that the lumped state space is the image of the lumping function and that all network states $x\in\mathcal{X}$ which are mapped to the same $y\in\mathcal{Y}$ will be aggregated.

Later in this section, we will discuss concrete realizations of $\mathcal{L}$. In particular, we will construct $\mathcal{L}$ based on a node partitioning and a counting abstraction of our choice. Next, we define the the transition rates $q(y,y^{\prime})$ (where $y,y^{\prime}\in\mathcal{Y}$, $y\neq y^{\prime}$) between the states of the lumped Markov chain:

[TABLE]

This is simply the mean transition rate at which an original state from $x$ goes to some $x^{\prime}\in\mathcal{L}^{-1}(y^{\prime})$. Technically, Eq. (1) corresponds to the following lumping assumption: we assume that at each point in time all network states belonging to a lumped state $y$ are equally likely.

4.1 Partition-Based Lumping

Next, we construct the lumping function $\mathcal{L}$. Because we want to make our lumping aware of the contact network’s topology, we assume a given partitioning $\mathcal{P}$ over the nodes $\mathcal{N}$ of the contact network. That is, $\mathcal{P}\subset 2^{\mathcal{N}}$ and $\bigcup_{P\in\mathcal{P}}P=\mathcal{N}$ and all $P\in\mathcal{P}$ are disjoint and non-empty. Based on the node partitioning, we can now impose different kinds of counting abstractions on the network state. This work considers two types: counting nodes and counting edges. The counting abstractions are visualized in Fig. 3. A full example of how a lumped CTMC of an $\mathtt{SI}$ model is constructed using the node-based counting abstraction is given in Fig. 2.

4.1.1 Node-Based Counting Abstraction

We count the number of nodes in each state and partition. Thus, for a given network state $x\in\mathcal{X}$, we use $y(s,P)$ to denote the number of nodes in state $s\in\mathcal{S}$ in partition $P\in\mathcal{P}$. The lumping function $\mathcal{L}$ projects $x$ to the corresponding counting abstraction. Formally:

[TABLE]

4.1.2 Edge-Based Counting Abstraction

Again, we assume that a network state $x$ and a node partitioning $\mathcal{P}$ are given. Now we count the edges, that is for each pair of states $s,s^{\prime}\in\mathcal{S}$ and each pair of partitions $P,P^{\prime}\in\mathcal{P}$, we count $y(s,P,s^{\prime},P^{\prime})$ which is the number of edges $(n,n^{\prime})\in\mathcal{E}$ where $x(n)=s$, $n\in P$, $x(n^{\prime})=s^{\prime}$, $n^{\prime}\in P^{\prime}$. Note that this includes cases where $P=P^{\prime}$ and $s=s^{\prime}$. However, only counting the edges does not determine how many nodes there are in each state (see Fig. 3 for an example).

In order to still have this information encoded in each lumped state, we slightly modify the network structure by adding a new dummy node $n_{\star}$ and connecting each node to it . The dummy node has a dummy state denoted by $\star$ which never changes, and it can be assigned to a new dummy partition $P_{\star}$. Formally,

[TABLE]

Note that the rate function $f$ ignores the dummy node. The lumped representation is then given as:

[TABLE]

4.1.3 Example

Fig. 2 illustrates how a given partitioning and the node-based counting approach induces a lumped CTMC. The partitions induced by the edge-based counting abstracting are also shown. In this example, the edge-based lumping aggregates only isomorphic network states.

4.2 Graph Partitioning

Broadly speaking, we have three options to partition the nodes based on local features (e.g., its degree) or global features (e.g., communities in the graph) or randomly. As a baseline, we use a random node partitioning. Therefore, we fix the number of partitions and randomly assign each node to a partition while enforcing that all partitions have, as far as possible, the same number of elements.

Moreover, we investigate a degree-based partitioning, where we define the distance between to nodes $n,n^{\prime}$ as their relative degree difference (similar to [28]):

[TABLE]

We can then use any reasonable clustering algorithm and build partitions (i.e., clusters) with the distance function. In this work, we focus on bottom-up hierarchical clustering as it provides the most principled way of precisely controlling the number of partitions. Note that, for the sake of simplicity (in particular, to avoid infinite distances), we only consider contact networks where each node is reachable from every other node. We break ties arbitrarily.

To get a clustering considering global features we use a spectral embedding of the contract network. Specifically, we use the spectral_layout function from the NetworkX Python-package [21] with three dimensions and perform hierarchical clustering on the embedding. In future research, it would be interesting to compute node distances based on more sophisticated graph embedding as the ones proposed in [16]. Note that in the border cases $|\mathcal{P}|=1$ and $|\mathcal{P}|=|\mathcal{N}|$ all methods yield the same partitioning.

5 Markov Population Models

Markov Population Models (MPMs) are a special form of CTMCs where each CTMC state is a population vector over a set of species. We use $\mathcal{Z}$ to denote the finite set of species (again, with an implicit enumeration) and $\mathbf{y}\in\mathbb{Z}_{\geq 0}^{|\mathcal{Z}|}$ to denote the population vector. Hence, $\mathbf{y}[z]$ identifies the number of entities of species $z$. The stochastic dynamics of MPMs is typically expressed as a set of reactions $\mathcal{R}$, each reaction, $(\alpha,\mathbf{b})\in\mathcal{R}$, is comprised of a propensity function $\alpha:\mathbb{Z}_{\geq 0}^{|\mathcal{Z}|}\rightarrow\mathbb{R}_{\geq 0}$ and a change vector $\mathbf{b}\in\mathbb{Z}^{|\mathcal{Z}|}$. When reaction $(\alpha,\mathbf{b})$ is applied, the system moves from state $\mathbf{y}$ to state $\mathbf{y}+\mathbf{b}$. The corresponding rate is given by the propensity function. Therefore, we can rewrite the transition matrix of the CTMC as333Without loss of generality, we assume that different reactions have different change vectors. If this is not the case, we can merge reactions with the same update by summing their corresponding rate functions.:

[TABLE]

Next, we show that our counting abstractions have a natural interpretation as MPMs.

5.1 Node-Based Abstraction

First, we define the set of species $\mathcal{Z}$. Conceptually, species are node states which are aware of their partition:

[TABLE]

Again, we assume an implicit enumeration of $\mathcal{Z}$. We use $z.s$ and $z.P$ to denote the components of a give species $z$.

We can now represent the lumped CTMC state as a single population vector $\mathbf{y}\in\mathbb{Z}_{\geq 0}^{|\mathcal{Z}|}$, where $\mathbf{y}[z]$ the number of nodes belonging to species $z$ (i.e., which are in state $z.s$ and partition $z.P$). The image of the lumping function $\mathcal{L}$, i.e. the lumped state space $\mathcal{Y}$, is now a subset of non-negative integer vectors: $\mathcal{Y}\subset\mathbb{Z}_{\geq 0}^{|\mathcal{Z}|}$.

Next, we express the dynamics by a set of reactions. For each rule $r=s_{1}\xrightarrow{f}s_{2}$ and each partition $P\in\mathcal{P}$, we define a reaction $(\alpha_{r,P},\mathbf{b}_{r,P})$ with propensity function as:

[TABLE]

where $\mathbf{m}_{x,n}$ denotes the neighborhood vector of $n$ in network state $x$. Note that this is just the instantiation of Equation 1 to the MPM framework.

The change vector $\mathbf{b}_{r,P}\in\mathbb{Z}^{|\mathcal{Z}|}$ is defined element-wise as:

[TABLE]

Note that $s_{1},s_{2}$ refer to the current rule and $z.s$ to the entry of $\mathbf{b}_{r,P}$.

5.2 Edge-Based Counting Abstraction

We start by defining a species neighborhood. The species neighborhood of a node $n$ is a vector $\mathbf{v}\in\mathbb{Z}_{\geq 0}^{|\mathcal{Z}|}$, where $\mathbf{v}[z]$ denotes the number of neighbors of species $z$. We define $\mathcal{V}_{n}$ to be the set of possible species neighborhoods for a node $n$, given a fixed contact network and partitioning. Note that we still assume that a dummy node is used to encode the number of states in each partition.

Assuming an arbitrary ordering of pairs of states and partitions, we define

[TABLE]

Let us define $\mathcal{V_{P}}$ to be the set of partition neighborhoods all nodes in $P$ can have:

[TABLE]

For each rule $r=s_{1}\xrightarrow{f}s_{2}$, and each partition $P\in\mathcal{P}$, and each $\mathbf{v}\in\mathcal{V}_{P}$, we define a propensity function $\alpha_{r,P,\mathbf{v}}$ with:

[TABLE]

Note that the propensity does not actually depend on $\mathbf{v}$, it is simply individually defined for each $\mathbf{v}$. The reason for this is that the change vector depends on the a node’s species neighborhood. To see this, consider a species $z=(s_{source},P_{source},s_{target},P_{target})$, corresponding to edges connecting a node in state $s_{source}$ and partition $P_{source}$ to a node in state $s_{target}$ and partition $P_{target}$. There are two scenarios in which the corresponding counting variable has to change: (a) when the node changing state due to an application of rule $r$ is the source node, and (b) when it is the target node. Consider case (a); we need to know how many edges are connecting the updated node (which was in state $s_{1}$ and partition $P$) to a node in state $s_{target}$ and partition $P_{target}$. This information is stored in the vector $\mathbf{v}$, specifically in position $\mathbf{v}[s_{target},P_{target}]$. The case in which the updated node is the target one is treated symmetrically. This gives rise to the following definition:

[TABLE]

The first two lines of the definition handle cases in which the node changing state is the source node, while the following two lines deal with the case in which the node changing state appears as target.

Fig. 4 illustrates how a lumped network state is influenced by the application of an infection rule.

5.3 Direct Construction of the MPM

Approximating the solution of an $\mathtt{SIS}$-type process on a contact network by lumping the CTMC first, already reduces the computational costs by many orders of magnitude. However, this scheme is still only applicable when it is possible to construct the full CTMC in the first place. Recall that the number of network states is exponential in the number of nodes of the contact network, that is, $|\mathcal{X}|=|\mathcal{S}|^{|\mathcal{N}|}$.

However, in recent years, substantial effort was dedicated to the analysis of very small networks [46, 23, 30, 33, 35]. One reason is that when the size of a network increases, the (macro-scale) dynamics becomes more deterministic because stochastic effects tend to cancel out. For small contact networks, however, methods which capture the full stochastic dynamics of the system, and not only the mean behavior, are of particular importance.

A substantial advantage of the reduction to MPM is the possibility of constructing the lumped CTMC without building the full CTMC first. In particular, this can be done exactly for the node counting abstraction. On the other hand, for the edge counting we need to introduce an extra approximation in the definition of the rate function, roughly speaking introducing an approximate probability distribution over neighboring vectors, as knowing how many nodes have a specific neighboring vector requires us full knowledge of the original CTMC. We present full details of such direct construction in Appendix 8.

5.4 Complexity of the MPM

The size of the lumped MPM is critical for our method, as it determines which solution techniques are computationally tractable and provides guidelines on how many partitions to choose. There are two notions of size to consider: (a) the number of population variables and (b) the number of states of the underlying CTMC. While the latter governs the applicability of numerical solutions for CTMCs, the former controls the complexity of a large number of approximate techniques for MPMs, like mean field or moment closure.

Node-based abstraction.

In this abstraction, the population vector is of length $|\mathcal{S}|\cdot|\mathcal{P}|$, i.e. there is a variable for each node state and each partition.

Note that the sum of the population variables for each partition $P$ is $|P|$, the number of nodes in the partition. This allows us to count easily the number of states of the CTMC of the population model: for each partition, we need to subdivide $|P|$ different nodes into $|\mathcal{S}|$ different classes, which can be done in $\binom{|P|+|\mathcal{S}|-1}{|\mathcal{S}|-1}$ways, giving a number of CTMC states exponential in the number $|\mathcal{S}|$ of node states and $|\mathcal{P}|$ of partitions, but polynomial in the number of nodes:

[TABLE]

Edge-based abstraction.

The number of population variables, in this case, is one for each edge connecting two different partitions, plus those counting the number of nodes in each partition and each node state, due to the presence of the dummy state. In total, we have $\frac{q(q-1)}{2}+q$ population variables, with $q=|\mathcal{S}|\cdot|\mathcal{P}|$.

In order to count the number of states of the CTMC in this abstraction, we start by observing that the sum of all variables for a given pair of partitions $P^{\prime},P^{\prime\prime}$ is the number of edges connecting such partitions in the graph. We use $\epsilon(P^{\prime},P^{\prime\prime})$ to denote the number of edges between $P^{\prime},P^{\prime\prime}$ (resp. the number of edges inside $P^{\prime}$ if $P^{\prime}=P^{\prime\prime}$). Thus,

[TABLE]

This is an over-approximation, because not all combinations are consistent with the graph topology. For example, a high number of infected nodes in a partition might not be consistent with a small number of $\texttt{I}-\texttt{I}$-edges inside the partition. Note that also this upper bound is exponential in $|\mathcal{S}|$ and $|\mathcal{P}|$ but still polynomial in the number of nodes $N$, differently from the original network model, whose state space is exponential in $N$.

The exponential dependency on the number of species (i.e., dimensions of the population vector) makes the explicit construction of the lumped state space viable only for very small networks with a small number of node states. However, this is typically the case for spreading models like SIS or SIR. Yet, also the number of partitions has to be kept small, particularly in realistic models. We expect that the partitioning is especially useful for networks showing a small number of large-scale homogeneous structures, as happens in many real-world networks [11].

An alternative strategy for analysis is to derive mean-field [4] or moment closure equations [39] for MPMs, which can be done without explicitly constructing the lumped (and the original) state space. These are sets of ordinary differential equation (ODE) describing the evolution of (moments of) the population variables. We refer the reader to [9] for a similar approach regarding the node-based abstraction.

6 Numerical Results

In this section, we compare the numerical solution of the original model—referred to as baseline model—with different lumped MPMs. The goal of this comparison is to provide evidence supporting the claim that the lumping preserves the dynamics of the original system, with an accuracy increasing with the resolution of the MPM. We will perform the comparison by solving numerically the ground and the lumped system, thus comparing the the probability of each state in each point in time. In practical applications of our method, exact transient or steady state solutions may not be feasible, but in this case we can still rely to approximation methods for MPM [4, 39]. Determining which of those techniques performs best in this context is a direction of future exploration.

A limit of the comparison based on numerical solution of the CTMC is that the state space of the original model has $|\mathcal{S}|^{|\mathcal{N}|}$ states, which limits the size of the contact network strongly444Code is available at github.com/gerritgr/Reducing-Spreading-Processes.

Let $P(X(t)=x)$ denote the probability that the baseline CTMC occupies network state $x\in\mathcal{X}$ at time $t\geq 0$. Furthermore, let $P(Y(t)=y)$ for $t\geq 0$ and $y\in\mathcal{Y}$ denote the same probability for a lumped MPM (corresponding to a specific partitioning and counting abstraction). To measure their difference, we first approximate the probability distribution of the original model using the lumped solution, invoking the lumping assumption which states that all network states which are lumped together have the same probability mass. We use $P_{L}$ to denote the lifted probability distribution over the original state space given a lumped solution. Formally,

[TABLE]

We measure the difference between the baseline and a lumped solution at a specific time point by summing up the difference in probability mass of each state, then take the maximum error in time:

[TABLE]

In our experiments, we used a small toy network with 13 nodes and 2 states ($2^{13}=8192$ network states). We generated a synthetic contact network following the Erdős–Rényi graph model with a connection probability of $0.5$. We use a SIS model with an infection rate of $\lambda=1.0$ and a recovery rate of $\mu=1.3$. Initially, we assign an equal amount of probability mass to all network states.

Fig. 5 shows the relationship between the error of the lumped MPM, the type of counting abstraction and the method used for node partitioning. We also report the mean difference together with the maximal difference over time.

From our results, we conclude that the edge-based counting abstraction yields a significantly better trade-off between state space size and accuracy. However, it generates larger MPM models than the node-based abstraction when adding a new partition. We also find that spectral and degree-based partitioning yield similar results for the same number of CTMC states and that random partitioning performed noticeably worse, for both edge-based and node-based counting abstractions.

7 Conclusions and Future Work

This work developed first steps in a unification of the analysis of stochastic spreading processes on networks and Markov population models. Since the so obtained MPM can become very large in terms of species, it is important to be able to control the trade-off between state space size and accuracy.

However, there are still many open research problems ahead. Most evidently, it remains to be determined which of the many techniques developed for the analysis of MPMs (e.g. linear noise, moment closure) work best on our proposed epidemic-type MPMs and how they scale with increasing size of the contact network. We expect also that these reduction methods can provide a good starting point for deriving advanced mean-field equations, similar to ones in [9]. Moreover, literature is very rich in proposed moment-closure-based approximation techniques for MPMs, which can now be utilized [42, 18]. We also plan to investigate the relationship between lumped mean-field equations [19, 28] and coarse-grained counting abstractions further.

Future work can additionally explore counting abstraction of different types, for instance, a neighborhood-based abstraction like the one proposed by James P. Gleeson in [12, 13].

Finally, we expect that there are many more possibilities of partitioning the contact network that remain to be investigated and which might have a significant impact on the final accuracy of the abstraction.

7.0.1 Acknowledgements

This research has been partially funded by the German Research Council (DFG) as part of the Collaborative Research Center \sayMethods and Tools for Understanding and Controlling Privacy. We thank Verena Wolf for helpful discussions and provision of expertise.

Appendix

8 Direct Construction of MPMs

Here, we prosper a way of directly deriving the lumped MPMs from the contact network without building the original CTMC first. We start with the node-based counting abstraction.

8.1 Node-Based Abstraction with General Rate Functions

Our general strategy is to iterate over the nodes in the contact network and to compute the mean rate attributed to that node over all $x\in\mathcal{X}$. Therefore, we consider the possible states of each node together with all possible species neighborhoods. The probability of a node $n$ being in state $s$ and having species neighborhood $\mathbf{v}$ is denoted as \Pr\big{(}X(n)=s,V(n)=\mathbf{v}\big{)}.

For a specific rule $r=s_{1}\xrightarrow{f}s_{2}$ and partition $P$, we can then describe $\alpha_{r,P}$ as:

[TABLE]

where $\mathbf{m}_{\mathbf{v}}$ is the neighborhood vector induced by $\mathbf{v}$, which we receive by grouping all partitions together. Note that it is not computationally necessary to actually iterate over all nodes in the partition. Instead we can group all nodes with the same partition neighborhood together, that is, all nodes $n^{\prime},n^{\prime\prime}\in P$ with $V_{n^{\prime}}=V_{n^{\prime\prime}}$ as the probability only depends on $\mathbf{v}$.

Computing the probability is the interesting part, we start by establishing that

[TABLE]

The first term in the product can be described by simply dividing the number of $s_{1}$-nodes in $P$ with the total number of nodes in $P$.

[TABLE]

The latter probability can be computed for each partition independently. This is because we know the number of nodes in each state in each partition. We also know that in partition $P$ the current node $n$ is already in state $s_{1}$, which we have to take into account. First, we define $\mathbf{y}_{P}\in\mathbb{Z}_{\geq 0}^{|\mathcal{S}|}$ to the the projection from $y$ to $P$. Thus, each entry is defined by:

[TABLE]

Likewise, we define $\mathbf{v}_{P}\in\mathbb{Z}_{\geq 0}^{|\mathcal{S}|}$, such that $\mathbf{v}_{P}[s]=\mathbf{v}[s,P]$. We also define $V(n)_{P}\in\mathbb{Z}_{\geq 0}^{|\mathcal{S}|}$ to be the number of neighbors of node $n$ in partition $P$ for each state. Finally, we define $\mathbf{y}_{P}^{s_{1}^{-}}$ to be the same vector as $\mathbf{y}_{P}$ except that the entry corresponding to state $s_{1}$ is subtracted by one (and truncated at zero). We can now rewrite the probability as:

[TABLE]

We use $p_{h}(\mathbf{k};\mathbf{K})$ to denote the probability mass function of the the multivariate hypergeometric distribution, where $\mathbf{k},\mathbf{K}$ denote to vectors over non-negative integers of the same length. That is, if $\mathbf{K}$ denotes the number of nodes in each state in a partition (resp., the number of marbles in an urn with different colors), then, $p_{h}(\mathbf{k};\mathbf{K})$ denotes the probability of drawing exactly $\mathbf{k}[s]$ nodes (resp. marbles) of each state (resp. color).

8.2 Reaction Networks and Linear Models

A special case of MPMs are biochemical reaction networks, where the species represent different types of molecules. The change vectors and corresponding propensity functions can elegantly be expressed as monomolecular ($\texttt{A}\rightarrow\texttt{B}$) and bimolecular ($\texttt{A}+\texttt{B}\rightarrow\texttt{C}+\texttt{D}$) reaction rules ($\texttt{A},\texttt{B},\texttt{C},\texttt{D}\in\mathcal{Z}$).

8.2.1 Reduction to Biochemical Reaction Networks

Most classical models in computational epidemiology are solely comprised of node-based rules (like the curing rule) and edge-based rules (like the infecting propagation rule). We call these linear models. Node-based rules, also referred to as spontaneous or independent rules, have a constant rate function, i.e., $f(\mathbf{m})=\mu$. Edge-based rules, also referred to as contact rules, are linear in exactly one dimension, i.e., they have the form $f(\mathbf{m})=\lambda\mathbf{m}[s]$.

Linear models are special because not the whole neighborhood is important for the rate of a rule but only the expected number of neighbors in a certain state. This makes the rule very similar to monomolecular and bimolecular reaction rates in MPMs. In fact, we can model the whole dynamics as a set of reaction over the species $\mathcal{Z}$.

Chemical reaction networks are a special case of Markov population models. In a chemical reaction network the state space is given by population vectors over species and molecular reactions have the form $\texttt{A}\xrightarrow{a}\texttt{C}$ or $\texttt{A}+\texttt{B}\xrightarrow{b}\texttt{C}+\texttt{D}$, where A,B,C,D denote species and $a,b\in\mathbb{R}_{\geq 0}$ are reaction rate constants.

For each node-based rule $s_{1}\xrightarrow{\mu}s_{2}$, we construct the reactions

[TABLE]

For each edge-based rule $s_{1}\xrightarrow{f}s_{1}$, $f(\mathbf{m})=\lambda\mathbf{m}[s^{\prime}]$, we construct the reactions

[TABLE]

where $w_{P,P^{\prime}}$ denotes the mean number of edges of a random node in $P$ with nodes in $P^{\prime}$, that is555Note that, despite the duple notation, we only count edges once:

[TABLE]

with

[TABLE]

8.3 Edge-Based Counting Abstraction

For each rule $r=s_{1}\xrightarrow{f}s_{2}$, and each partition $P\in\mathcal{P}$, and each $\mathbf{v}\in\mathcal{V}_{P}$, we define a propensity function $\alpha_{r,P,\mathbf{v}}$ with:

[TABLE]

Again, we use

[TABLE]

to compute this probability, where we can solve \Pr\big{(}X(n)=s_{1}\big{)} exactly as before.

Since we have now information about the edges, we can derive the probability of neighborhoods more precisely. In fact, we can directly construct the set of candidate neighbors from $\mathbf{y}$. Therefore, we define a vector $\mathbf{y}_{s,P,P^{\prime}}\in\mathbb{Z}_{\geq 0}^{|\mathcal{S}|}$, where entry $\mathbf{y}_{s,P,P^{\prime}}[s^{\prime}]$ specifies the number of neighbors of a random node in state $s$ and partition $P$, which lie in partition $P^{\prime}$ and occupy state $s^{\prime}$. Formally:

[TABLE]

This gives rise to the final approximation of the probability of neighborhood species:

[TABLE]