Exact analysis of summary statistics for continuous-time discrete-state Markov processes on networks using graph-automorphism lumping

We propose a unified framework to represent a wide range of continuous-time discrete-state Markov processes on networks, and show how many network dynamics models in the literature can be represented in this unified framework. We show how a particular sub-set of these models, referred to here as single-vertex-transition (SVT) processes, lead to the analysis of quasi-birth-and-death (QBD) processes in the theory of continuous-time Markov chains. We illustrate how to analyse a number of summary statistics for these processes, such as absorption probabilities and first-passage times. We extend the graph-automorphism lumping approach [Kiss, Miller, Simon, Mathematics of Epidemics on Networks, 2017; Simon, Taylor, Kiss, J. Math. Bio. 62(4), 2011], by providing a matrix-oriented representation of this technique, and show how it can be applied to a very wide range of dynamical processes on networks. This approach can be used not only to solve the master equation of the system, but also to analyse the summary statistics of interest. We also show the interplay between the graph-automorphism lumping approach and the QBD structures when dealing with SVT processes. Finally, we illustrate our theoretical results with examples from the areas of opinion dynamics and mathematical epidemiology.

Dynamical processes on networks are one of the main topics in network science (Barrat et al. 2008; Castellano et al. 2009; Newman 2003; 2010; Porter and Gleeson 2016). Numerous applications have been modelled as dynamical processes on networks, including epidemics (Kiss et al. 2017; Pastor-Satorras et al. 2015), magnetism (Glauber 1963), opinion dynamics (Galam 2002; Sood and Redner 2005; Sznajd-Weron and Sznajd 2000), diffusion of innovations (Bass 1969; Mellor et al. 2015; Melnik et al. 2013; Watts 2002), rumour spread (Daley and Kendall 1964; Goldenberg et al. 2001; Kempe et al. 2003), meme popularity (Gleeson et al. 2014), cultural polarisation (Axelrod 1997; Castellano et al. 2000), racial segregation (Schelling 1969; 1971), stock market trading (Kirman 1993), cascading failures (Gleeson et al. 2012; Haldane and May 2011; Motter and Lai 2002) and language evolution (Baronchelli et al. 2006; Bonabeau et al. 1995; Castelló et al. 2006). The mathematical analysis of such models has highlighted the rich yet subtle dependence of dynamical phenomena on network topology (Newman 2010; Porter and Gleeson 2016; Durrett 2007). However, the standard mathematical approach is via mean-field theories (Kiss et al. 2017), where dynamical correlations, clustering and modularity are ignored (Gleeson et al. 2012). While it is possible to derive high-accuracy mean-field approximations (Gleeson 2011) and some approximation approaches are amenable to rigorous analysis (Li et al. 2012; Van Mieghem 2011), there is still no theory that quantifies the accuracy of mean-field approximations or characterises the types of networks or dynamical processes where approximations work well.

In contrast to the use of approximate mean-field theories or simulation studies (Stoll et al. 2012), recent approaches in mathematical epidemiology have focused on the exact analysis of infection spread dynamics occurring on small networks, for example to quantify the importance of nodes, in terms of outbreak size, vaccination and early infection in SIR epidemics (Holme 2017), and to compute SIS extinction times using computational algebra for all sufficiently small graphs (Holme and Tupikina 2018). By focusing on small networks, it is possible to include heterogeneous rates of infection and recovery in the context of particular applications, such as the spread of hospital-acquired infections in intensive care units (López-García 2016), and to analyse these systems in terms of a number of performance measures. Thus, the aim is usually to compute summary statistics related to the dynamical process (Economou et al. 2015), instead of focusing on analysing the complete transient dynamics of the process, which are usually more complex to study (Keeling and Ross 2009).

Loosely speaking, a dynamical process on a network describes how the state of each node or vertex in the network changes in time. In this paper, we focus on dynamical processes that are described in terms of continuous-time Markov chains (CTMCs), although other types of dynamical processes are possible (for example where the state-space is continuous (Deffuant et al. 2000; Fortunato 2004; Hegselmann and Krause 2002), and/or the dynamics deterministic (Nakao and Mikhailov 2010; Rodrigues et al. 2016; Ward and Grindrod 2014) or non-Markovian (Kiss et al. 2015; Castro et al. 2018)). Thus for brevity, we will refer here to dynamical processes on networks described by CTMCs as simply network dynamics. It is well-known that the state-space of network dynamics grows exponentially with the number of nodes, so can be extremely large for even fairly small networks. Furthermore, network dynamics are determined by the infinitesimal generator matrix, which contains the rates at which the process moves between pairs of states. Given the large size of this matrix, even constructing it can prove challenging (Van Mieghem et al. 2009). In the area of mathematical epidemiology it has been recently shown how many network dynamical processes can be represented as quasi-birth-and-death processes (QBDs) (Simon et al. 2011; López-García 2016; Economou et al. 2015) by appropriately ordering the states within the state-space of the corresponding CTMC, so that its infinitesimal generator matrix is significantly sparse and tridiagonal-by-blocks. QBDs exploit the tri-diagonal structure of transitions occurring across groups or levels of states, and summary statistics in these systems can be analysed by means of computationally efficient algorithms (Latouche and Ramaswami 1999). While this approach is standard in areas such as queuing theory (Bean et al. 1997; Gómez-Corral and López-García 2014) and population dynamics (Gómez-Corral and López-García 2015), there is potential for it to be adopted more widely in applications of dynamics on networks in other areas of network science.

One way to reduce the size of the state-space is via lumping (Kemeny and Snell 1960), where several states in the state-space are combined together into one lumped state. Under certain conditions, quantities of interest about the unlumped system, such as its stationary distribution or transient probabilities, can be determined from the lumped system (Buchholz, P 1994). Efficient algorithms have been developed that can construct the optimal lumping quotient in \(\mathcal {O}(m\log (n))\) time, where n is the number of states and m is the number of transitions (Valmari and Franceschinis 2010; Derisavi et al. 2003). However, for network dynamics the number of states is exponential in the number of vertices N, and the number of transitions is typically of order N times the number of states. It has also been recently shown that network symmetries can be used to lump network dynamics (Kiss et al. 2017; Simon et al. 2011). For studies of small networks, such as those described above, lumping state-space on the basis of even relatively few symmetries can dramatically increase the size of the network that can be studied.

The main contribution of this paper is to illustrate how lumping based on network symmetries can be expressed in a matrix-oriented fashion for a broad class of network dynamics. This includes a large sub-class that give rise to QBD structures for which there are efficient algorithmic and matrix-oriented methods for the exact analysis of summary statistics. In the “CTMC dynamics on networks” section we present a unified CTMC framework that captures a wide range of network dynamics, the typical summary statistics of interest and how network symmetries can be used to lump the state-space using a matrix-oriented approach. In the “Single-vertex transition processes” section we focus on a particular sub-class of models that we call single-vertex transition processes that result in QBDs, and illustrate how the summary statistics described in the “CTMC dynamics on networks” section can be computed efficiently both for the full and lumped state-spaces. In the “Applications to classes of model” section, we show how our theoretical approach applies to various types of models using specific examples from opinion dynamics and mathematical epidemiology.

We start by introducing the general set-up of the type of CTMC network dynamics that we consider. This sets out our notation and crucially the structure of state-space. Model specific transition rates are discussed in detail in “Applications to classes of model” section.

Let us consider a network \({\mathcal {G}}=({\mathcal {V}},{\mathcal {E}})\) where \({\mathcal {V}}=\{1,\dots,N\}\) is the set of vertices and \({\cal E}\subseteq {\mathcal {V}}\times {\mathcal {V}}\) is the set of edges of the network. In this paper we focus on simple networks, i.e. unweighted, undirected networks with no self- or multi-edges. We consider that, over time, vertices can take a number of different vertex-states from the set \({\mathcal {W}}=\{1,\dots,M\}\), and that vertices change their vertex-state according to a continuous-time Markov chain (CTMC) \({\mathcal {X}}=\{\textbf {X}(t):\ t\geq 0\}\). The state of \({\mathcal {X}}\) at time t≥0, X(t), is a random vector \(\textbf {X}(t)=\left (X_{1}(t),\dots,X_{N}(t)\right)\) where \(X_{i}(t)\in \mathcal {W}\) represents the vertex-state of vertex i at time t≥0, for 1≤i≤N. Since each vertex takes vertex-states in \({\mathcal {W}}\), it is clear that the space of states of the CTMC is given by \({\mathcal {S}}={\cal W}^{\mathcal {V}}=\{\left (w_{1},\dots,w_{N}\right):\ w_{j}\in {\mathcal {W}},\ 1\leq j\leq N\}\), with cardinality \(\#{\mathcal {S}}=M^{N}\), so that it can be written as \({\cal S}=\{S_{1},S_{2},\dots,S_{M^{N}}\}\), with any arbitrary order of states. That is, each state \(S\in {\mathcal {S}}\) of the process \({\mathcal {X}}\) is just a permutation with repetition of N vertex-states.

For a given state \(S\in {\mathcal {S}}\), we denote by \(S(v)\in {\mathcal {W}}\) the vertex-state of vertex v in state S, for any 1≤v≤N. Also, given two states \(S_{i},S_{j}\in {\mathcal {S}}\), the sub-set \(U(S_{i},S_{j})\subset {\mathcal {V}}\) of so-called transition vertices contains those vertices u∈U(S_i,S_j) such that S_i(u)≠S_j(u) and S_i(v)=S_j(v) for every \(v\in {\mathcal {V}}\setminus U(S_{i},S_{j})\) (that is, those vertices that differ in their vertex-state between states S_i and S_j). For the (one-jump) transition between states \(S_{i}\rightarrow S_{j}\), we denote the transition rate by \(q_{S_{i}\rightarrow S_{j}}\geq 0\). We will discuss how the transition rates \(q_{S_{i}\rightarrow S_{j}}\) can be determined for specific types of model in “Applications to classes of model” section. For any given state \(S_{i}\in {\mathcal {S}}\), we denote by \(A(S_{i})\subset {\mathcal {S}}\) the set of states directly accessible (that is, in one jump of the process) from S_i; thus \(A(S_{i})=\{S_{j}\in {\mathcal {S}}:\ q_{S_{i}\rightarrow S_{j}}>0\}\).

The dynamics of the CTMC \({\mathcal {X}}\) can be then analysed in terms of its transition probabilities \(p_{S_{j},S_{i}}(t)=\mathbb {P}(\textbf {X}(t)=S_{i}\ |\ \textbf {X}(0)=S_{j})\), \(\forall S_{i},S_{j}\in {\mathcal {S}}\) and t≥0. These probabilities, which can be stored in a matrix \(\textbf {P}(t)=(p_{S_{j},S_{i}}(t))_{i\in \{1,\dots,M^{N}\},\ j\in \{1,\dots,M^{N}\}}\), can be computed by solving the master equation

where Q is the infinitesimal generator of \({\mathcal {X}}\). In particular, the matrix Q contains the transition rates between states in an ordered fashion, as follows

where elements in the diagonal of Q are just the negative sum of elements in each row.

For an arbitrary \(m\in \mathcal {W}\), say here \(m\in \mathcal {W}\) without any loss of generality, one can organise the space of states \({\mathcal {S}}\) into groups or levels,

That is, each level \({\mathcal {S}}^{(i)}\) is formed by all the states in \({\mathcal {S}}\) that contain exactly i vertices with vertex-state M. Each level can be further divided into sub-levels based on different vertex-states (Economou et al. 2015), although this will not be necessary in this paper. It is clear that the number of states in level \({\mathcal {S}}^{(i)}\) is the combinatorial number \(\#{\mathcal {S}}^{(i)}=\binom {N}{i}(M-1)^{N-i}=\frac {N!}{i!(N-i)!}(M-1)^{N-i}\), which is consistent since

One can denote then by \({\mathcal {S}}^{(i)}=\left \{S^{(i)}_{1},\dots,S^{(i)}_{\#{\mathcal {S}}^{(i)}}\right \}\), for 0≤i≤N, the states contained inside each level. Here, \(S^{(i)}_{j}\) denotes the j^th state inside level \({\mathcal {S}}^{(i)}\), according to some arbitrary order, for \(1\leq j\leq \#{\mathcal {S}}^{(i)}\). Then, by ordering states in \({\mathcal {S}}\) in terms of levels, \({\mathcal {S}}^{(0)}\prec {\mathcal {S}}^{(1)}\prec \dots \prec {\mathcal {S}}^{(N)}\), the matrix Q can be written as

where each sub-matrix Q_ij contains the infinitesimal transition rates from states in level \({\mathcal {S}}^{(i)}\) to states in level \({\mathcal {S}}^{(j)}\). We note here that \({\mathcal {S}}^{(N)}=\{(M,\dots,M)\}\) has a single state, so that the matrices Q_Nj are in fact row vectors, matrices Q_jN are column vectors, and Q_NN is a single-element matrix.

Instead of analysing the master Eq., (1), directly, one can analyse the dynamics of the process \({\mathcal {X}}\) by focusing on a number of summary statistics or stochastic descriptors. In the area of mathematical epidemiology, and when analysing the spread of a disease through a population by means of network dynamics, these summary statistics could represent the size of the outbreak (Economou et al. 2015), the time until the epidemic dies out (Economou et al. 2015), the area under the curve of infectives (Ball and Clancy 1995), or the probability of infection for a particular node of interest (López-García 2016). We refer the reader to (Keeling and Ross 2007) for a discussion on how exact and matrix-oriented approaches can be implemented to compute summary statistics of interest in continuous-time Markov chains.

Summary statistics

We focus in this section on quantities of interest related to absorption. That is, the arrival of the dynamical process to a particular absorbing state, so that the process remains indefinitely in this state after arrival. However, we point out that our arguments would be easily generalisable to other situations, such as the computation of steady-state probabilities for positive recurrent CTMCs (Buchholz, P 1994; Kulkarni 2016), or other summary statistics or quantities of interest (López-García 2016; Economou et al. 2015), depending on the particular network dynamics under study.

Probabilities of absorption

Let us consider that the process \({\mathcal {X}}\) has two absorbing states, so that

where S_a and S_b are two arbitrary absorbing states and \({\mathcal {C}}=\{S_{1},S_{2},\dots,S_{M^{N}-2}\}\) is the class of transient states. An example of this situation is provided in the “Biased voter model” section. We note that although we consider here two absorbing states in our system, our arguments generalise to processes with more than two absorbing states. We are interested in computing the probabilities of absorption into S_a from an initial state S_i,

Probabilities \(\{p_{S_{i}}^{(a)}:\ S_{i}\in {\mathcal {C}}\}\) can be computed by following a first-step argument, which yields the system of linear equations

where \(\Delta _{S_{i}}=\sum \limits _{S_{k}\in A(S_{i})}q_{S_{i}\rightarrow S_{k}}\), and with boundary conditions \(p^{(a)}_{S_{a}}=1\) and \(p^{(a)}_{S_{b}}=0\). Since each equation in (4) corresponds to an initial transient state \(S_{i}\in {\mathcal {C}}\), the system (4) can be written in matrix form as

For \(S_{i},S_{j}\in \mathcal {C}\), the ijth component of the matrix A is the one-step jump transition probability \(q_{S_{i}\rightarrow S_{j}}/\Delta _{S_{i}}\), the ith component of the column vector p^(a) is the probability \(p_{S_{i}}^{(a)}\) and the ith component of the column vector b^(a) is the one-step jump transition probability \(q_{S_{i}\rightarrow S_{a}}/\Delta _{S_{i}}\), which might be zero if S_a∉A(S_i).

Time until absorption

Let us now consider the alternative situation where the process \({\mathcal {X}}\) has only one (arbitrary) absorbing state, S_a, so that

where \({\mathcal {C}}=\{S_{1},S_{2},\dots,S_{M^{N}-1}\}\) is the class of transient states. We note that if one has several absorbing states in the process, and is interested in the time until absorption (into any absorbing state), one way to proceed is to combine all these absorbing states into an absorbing macro-state, and study the time until absorption into this macro-state.

Here, \(\lim \limits _{t\rightarrow +\infty }\mathbb {P}(\textbf {X}(t)=S_{a})=1\) regardless of the initial state X(0). Our interest is instead in analysing the time until absorption, \(T=\inf \{t\geq 0:\ \textbf {X}(t)=S_{a}\}\). To this end, we can define its Laplace-Stieltjes transform as \(\phi _{S_{i}}(z)=E[e^{-zT}\ |\ \textbf {X}(0)=S_{i}]\), for \(\Re (z)\geq 0\) and \(S_{i}\in {\mathcal {S}}\), and by a first-step argument one can obtain the system of linear equations

with boundary conditions \(\phi _{S_{a}}(z)=1\) for all \(\Re (z)\geq 0\). This system of equations allows one to compute any order moment of T by direct differentiation, since given any initial state \(S_{i}\in {\mathcal {S}}\), \(m^{(k)}_{S_{i}}=E[T^{k}\ |\ \textbf {X}(0)=S_{i}]=(-1)^{k}\left.\frac {d^{k}\phi _{S_{i}}(z)}{dz^{k}}\right |_{z=0}\). For example, for k=1 the first order moments \(m^{(1)}_{S_{i}}=E[T\ |\ \textbf {X}(0)=S_{i}]\) can be computed from the system of equations

with boundary conditions \(m^{(1)}_{S_{a}}=0\). The system of equations given by (7) can be expressed in matrix form as

where m⁽¹⁾ is a column vector containing the mean times \(m^{(1)}_{S_{i}}\) for the different initial states \(S_{i}\in {\mathcal {C}}\). For \(S_{i},S_{j}\in \mathcal {C}\), the ijth component of the matrix B is the one-step jump transition probability \(q_{S_{i}\rightarrow S_{j}}/\Delta _{S_{i}}\) and the ith component of the column vector c is \(1/\Delta _{S_{i}}\).

In fact, the time to absorption T is known to follow a phase-type distribution (He 2014, Chapter 1) PH(τ,T), where instead of considering a single initial state \(S\in {\mathcal {S}}\), one can consider a row vector τ of initial probabilities for states \(S\in {\mathcal {C}}\), and T is the sub-matrix of the infinitesimal generator corresponding to transient states in \({\mathcal {C}}\). The distribution function of T, for any t≥0, is given by (He 2014, Definition 1.2.1)

where 1 is a column vector of ones, and any k^th order moment of T can be obtained as (He 2014, Proposition 1.2.2)

We note that for k=1, this is equivalent to solving system (8).

The main problem when dealing with (1) is its cardinality, where the number of ODEs is in correspondence with the number of states, M^N. This same problem affects to the computation of the summary statistics of interest, and in particular systems (4) and (7). Building on recent work (Kiss et al. 2017; Simon et al. 2011; Ward and Evans 2019), we now describe how network symmetries can be used to reduce the size of the CTMC state-space and how this affects the summary statistics described in “Summary statistics” section.

A partition of state-space \({\mathcal {L}}=\{L_{1},L_{2},\dots,L_{r}\}\) is called a lumping if, roughly speaking, one can represent the original dynamics across states in \({\mathcal {S}}\), with dynamics across states in \({\mathcal {L}}\), while preserving the Markov property (Kemeny and Snell 1960). We define the collector matrix\(\textbf {V}\in \{0,1\}^{M^{N}\times r}\) (Buchholz, P 1994) whose ijth component is

More precisely, a necessary and sufficient condition for \({\mathcal {L}}\) to be a lumping (Kemeny and Snell 1960) is that for each L_i and L_j, there exists a \(\hat {q}_{L_{i}\rightarrow L_{j}}\) that satisfies

Also, following (Buchholz, P 1994), if \({\mathcal {L}}\) is a lumping then we say that it is a strict lumping if for each L_i and L_j, there also exists a \(\tilde {q}_{L_{i}\rightarrow L_{j}}\) that satisfies

We now show how one can derive a lumping of the state-space from symmetries of the network, by extending the arguments in (Kiss et al. 2017; Simon et al. 2011; Ward and Evans 2019), while providing a matrix-oriented representation of the lumping. An automorphism of a network \({\mathcal {G}}\) is a bijection \(g:{\mathcal {V}}\rightarrow {\mathcal {V}}\) such that \((u,v)\in {\mathcal {E}}\) if and only if \((g(u),g(v))\in {\mathcal {E}}\). We use the shorthand gv=g(v). The set of automorphisms of a network \({\mathcal {G}}\) form a permutation group \(H=\text {Aut}({\mathcal {G}})\) called the automorphism group of \({\mathcal {G}}\) (Godsil and Royle 2013). Automorphisms of the network permute vertices, so we need to define how they act on the state-space. Let H be the automorphism group of \({\mathcal {G}}\), then we define the action of g∈H on a state \(S_{i}\in \mathcal {S}\) to be

i.e. the vertex-state of v in gS_i is the same as the vertex-state of g⁻¹v in S_i. It is easy to prove that this is indeed an action of H on \(\mathcal {S}\) in the group theoretic sense. It follows from our definition (10) of the action of H on state-space that \({\mathcal {S}}\) is an H-set (Fraleigh 2003). This means that we can partition the state-space into equivalence classes where a pair of states \(S_{i},S_{j}\in {\mathcal {S}}\) are equivalent if and only if there is a g∈H such that S_i=gS_j. This partition is known as the orbit partition of the state-space.

Theorem 1.

Let \(\mathcal {G}=(\mathcal {V},\mathcal {E})\) be a finite network with automorphism group H and \(\mathcal {W}\) be a non-empty and finite set of vertex-states. Let \(\mathcal {S}=\mathcal {W}^{\mathcal {V}}\) be the state-space of a CTMC with infinitesimal generator Q, as defined in (2), with the action of H on \(\mathcal {S}\)defined in (10). If there is a subgroup G≤H such that

then the orbit partition of \(\mathcal {S}\) under the action of G is a strict lumping of \(\mathcal {S}\).

Proof

Let \({\mathcal {L}}=\{L_{1},L_{2},\dots,L_{r}\}\) be the orbit partition of the state-space \(\mathcal {S}\) under the action of G. Suppose that S_k∈L_i and S_l∈L_j. For any S_m∈L_i, we can find a g∈G such that S_m=gS_k; thus let S_n=gS_l∈L_j. By assumption, \(q_{S_{k}\rightarrow S_{l}}=q_{S_{m}\rightarrow S_{n}}\). Since \(\mathcal {S}\) is a G-set, it follows that g is a bijection on L_j. Thus for all S_k,S_m∈L_i we have

A similar argument can be used to show that for all S_l,S_n∈L_j we have

Consequently \({\mathcal {L}}\) is a strict lumping of \(\mathcal {S}\). □

Thus we can make use of network automorphism lumping whenever (11) is satisfied. It has been shown recently that a broad class of ‘homogeneous’ network dynamics models satisfy this criterion for the whole automorphism group H (Kiss et al. 2017; Simon et al. 2011; Ward and Evans 2019). In this paper, we go further and show that certain ‘heterogeneous’ models that respect (at least some) symmetries of the network also satisfy (11). It is straightforward to prove that a sub-set of H that satisfies (11) forms a group, and so this can be used to lump the state-space in accordance with Theorem 1. In what follows, we will assume that the models under consideration satisfy (11) (that is, that there are some symmetries in the network that can be exploited for lumping purposes).

We will now show how the lumping of the state-space reduces the dimensionality of the master equation (Kiss et al. 2017; Simon et al. 2011), and how this can be done in a matrix-oriented fashion. We define the distributor matrix\(\textbf {W}\in \mathbb {R}^{r\times M^{N}}\) (Buchholz, P 1994), whose ijth component is

It is easy to show that WV=I_r, where I_r is the r×r identity matrix, and

Furthermore, VW commutes with Q, which follows from the fact that \({\mathcal {L}}\) is a strict lumping of \({\mathcal {S}}\). To see this, note that for any \(L_{i},L_{j}\in \mathcal {L}\) we have

But for any S_k∈L_i and S_l∈L_j we have

thus it follows that QVW=VWQ.

Denote the lumped transition probability by \(\hat {p}_{L_{j},L_{i}}(t)=\mathbb {P}(\textbf {X}(t)\in L_{i}\ |\ \textbf {X}(0)\in L_{j})\), \(\forall L_{i},L_{j}\in {\mathcal {L}}\) and t≥0, and let \(\hat {\mathbf {P}}(t)= (\hat {p}_{L_{j},L_{i}}(t))_{i\in \{1,\dots,r\},\ j\in \{1,\dots,r\}}\). Thus \(\hat {\mathbf {P}}=\textbf {W}\textbf {P}\textbf {V}\) and with the properties described above it is easy to derive the lumped master equation

where \(\hat {\mathbf {Q}}=\mathbf {W}\mathbf {Q}\mathbf {V}\) is the lumped infinitesimal generator. Note that one must use the same ordering of states in \({\mathcal {S}}\) when constructing the matrices W, Q and V; see the example in the “Applications to classes of model” section.

To differentiate between the CTMCs on the full state-space and the lumped state-space, we will refer to the states of the lumped state-space as orbits in the remainder of this paper. This highlights the fact that the lumped state-space results from the orbit partition of the unlumped state-space under the action of a sub-group of the network automorphism group, as per Theorem 1.

Lumping and summary statistics

Lumping and probabilities of absorption

We now return to the situation described in the “Probabilities of absorption” section where we assume the state-space has two absorbing states S_a and S_b, and the set of transient states is \({\mathcal {C}}\). We assume that the subgroup G of the network \({\mathcal {G}}\) satisfying (11) is non-trivial, and that \({\mathcal {L}}=\{L_{1},L_{2},\dots,L_{r}\}\) is the orbit partition of \(\mathcal {C}\) under the action of G. We also assume that S_a and S_b are fixed by all elements in G so that they cannot be lumped together. It follows from the definition of lumping that for every \(L_{i}\in \mathcal {L}\) and S_k,S_l∈L_i, \(\Delta _{S_{k}}=\Delta _{S_{l}}\), thus for any S_k∈L_i we define \(\Delta _{L_{i}}=\Delta _{S_{k}}\).

Let \(\hat {p}_{L_{i}}^{(a)}\) be the probability of absorption into S_a from the lumping cell (i.e., orbit) \(L_{i}\in {\mathcal {L}}\). These probabilities can be obtained from the lumped transition rates \(\hat {q}_{L_{i}\rightarrow L_{j}}\) in the manner described in the “Probabilities of absorption” section, resulting in the system of equations

For \(L_{i},L_{j}\in \mathcal {L}\), the ijth component of the matrix \(\hat {\mathbf {A}}\) is the one-step jump transition probability \(\hat {q}_{L_{i}\rightarrow L_{j}}/\Delta _{L_{i}}\), the ith component of the column vector \(\hat {\mathbf {p}}^{(a)}\) is the probability \(\hat {p}_{L_{i}}^{(a)}\) and the ith component of the column vector \(\hat {\mathbf {b}}^{(a)}\) is the one-step jump transition probability \(\hat {q}_{L_{i}\rightarrow S_{a}}/\Delta _{L_{i}}\), which might be zero if S_a cannot be accessed in one jump from L_i.

We now show that the probabilities of absorption are the same for all states within the same lumping cell. An automorphism of the network g∈G is a permutation of the vertices in \(\mathcal {V}\), thus with the action (10) we can define Π_g to be the representation of g as a permutation matrix on the transient states \(\mathcal {C}\). When (11) is satisfied, we have Π_gA=AΠ_g and \(\mathbf {\Pi }_{g}\mathbf {b}^{(a)}=\textbf {b}^{(a)}\). Consequently

and hence \(\textbf {p}^{(a)}=\mathbf {\Pi }_{g}\textbf {p}^{(a)}\). Since the solution to \((\textbf {I}_{\#\mathcal {C}}-\textbf {A})\textbf {p}^{(a)}=\textbf {b}^{(a)}\) is unique, it follows that for any \(L_{i}\in \mathcal {L}\) and S_k,S_l∈L_i, \(\hat {p}_{L_{i}}^{(a)}=p_{S_{k}}^{(a)}=p_{S_{l}}^{(a)}\).

Lumping and time to absorption

We now consider the time to absorption for the situation described in the “Time until absorption” section, where there is one absorbing state S_a and the set of transient states is \({\mathcal {C}}\). If one orders states as

the infinitesimal generator is given by

Again we assume that the subgroup G of the network \({\mathcal {G}}\) is non-trivial, and that \({\mathcal {L}}=\{L_{1},L_{2},\dots,L_{r}\}\) is the orbit partition of \(\mathcal {C}\) under the action of G. Necessarily, S_a is fixed by all elements in G, and we can consider L_r=S_a is the orbit with this single absorbing state. Let V and W be the collector and distributor matrices of \(\mathcal {C}\) respectively, which have the form

Recall that the distribution of the time to absorption T is given by a phase-type distribution PH(τ,T), where τ is a row vector of the initial probabilities of being in states in \(\mathcal {C}\) and T is the matrix of transition rates between states in \(\mathcal {C}\), as described above. It is well known that a phase-type distribution is not uniquely defined by the pair (τ,T), and we will show that the distribution function of the time to absorption T is the same for the full and lumped Markov chains. Let \(\boldsymbol {\mu }=\boldsymbol {\tau }\tilde {\mathbf {V}}\) and \(\textbf {U}=\tilde {\mathbf {W}}\textbf {T}\tilde {\mathbf {V}},\) so U is the matrix of lumped one-jump transition rates between (transient) states in \(\{L_{1},L_{2},\dots,L_{r-1}\}\). Let 1_k be the column vector with k elements all equal to 1, then \(\mathbf {1}_{r}=\tilde {\mathbf {W}}\mathbf {1}_{\#\mathcal {C}}\) and \(\tilde {\mathbf {V}}\tilde {\mathbf {W}}\mathbf {1}_{\# \mathcal {C}}=\mathbf {1}_{\# \mathcal {C}}\). An additional consequence of the fact that VW commutes with Q, and hence T, is that \((\tilde {\mathbf {W}} \textbf {T}\tilde {\mathbf {V}})^{n}=\tilde {\mathbf {W}}\textbf {T}^{n}\tilde {\mathbf {V}}\), which is easily proved by induction. It follows that \(\exp (\textbf {U}t)=\tilde {\mathbf {W}}\exp (\textbf {T}t)\tilde {\mathbf {V}}\), \(\tilde {\mathbf {V}}\tilde {\mathbf {W}}\) commutes with \(\exp (\textbf {T}t)\) and consequently

Thus the distribution of time to absorption into S_a is the same for both (τ,T) and (μ,U). It follows that all moments of T will be the same for both the full and lumped state-spaces. Moreover the mean time to absorption for any pair of states in the same lumping cell are the same, which can be shown in a manner analogous to that used to show that the probabilities of absorption are the same.

We focus in this section on what we call single-vertex-transition (SVT) processes, which satisfy the condition that for every pair of states \(S_{i}\in {\mathcal {S}}\) and S_j∈A(S_i), U(S_i,S_j) has cardinality #U(S_i,S_j)=1. This means that transitions between states, according to the CTMC \({\mathcal {X}}\), can only occur between states S_i and S_j that differ in one vertex-state (that is, in each single jump, only one vertex-state is updated). If the states in \(\mathcal {S}\) are ordered by levels, as described in the “CTMC dynamics on networks” section, then it is clear that non-zero infinitesimal transition rates can only occur between any given state \(S\in {\mathcal {S}}^{(i)}\) and states in either the same level \({\mathcal {S}}^{(i)}\) (provided that \(\#\mathcal {W}>2\)), or states in adjacent levels\({\mathcal {S}}^{(i-1)}\) and \({\mathcal {S}}^{(i+1)}\). This means that (3) becomes

so that Q is tridiagonal-by-blocks, leading to a quasi-birth-and-death (QBD) process (Latouche and Ramaswami 1999).

Probabilities of absorption

Let us now consider the same situation as the one analysed in the “Probabilities of absorption” section, where the process \({\mathcal {X}}\) has two absorbing states. For convenience, let us assume that the two absorbing states are \(S_{a}=S^{(0)}_{1}\) (i.e., the first state in level \({\mathcal {S}}^{(0)}\)) and \(S_{b}=(M,\dots,M)\) (i.e., the only state in level \({\mathcal {S}}^{(N)}\)), so that

where \({\mathcal {C}}\) is the class of transient states. The motivation for assuming these particular choices of S_a and S_b is purely to place them at the beginning and end respectively of the order of states in state-space. Given the structure by levels above, it is clear that

where \({\cal \hat S}^{(0)}=\{S^{(0)}_{2},\dots,S^{(0)}_{\#{\mathcal {S}}^{(0)}}\}\) is just the level \({\mathcal {S}}^{(0)}\) without the absorbing state \(S^{(0)}_{1}\).

Since the state \(S^{(0)}_{1}\) is absorbing, \(A(S^{(0)}_{1})=\emptyset \), thus we note that

and Q_NN−1=0 and Q_NN=0, since state \((M,\dots,M)\) is also absorbing. The infinitesimal generator of the process \({\mathcal {X}}\) in this scenario is then

It is clear that the matrix \(\hat {\mathbf {Q}}\) contains the transition rates between states in the transient class \({\mathcal {C}}\), while the column vectors q^(a) and q^(b) contain the transition rates from states in the transient class \({\mathcal {C}}\) to the absorbing states \(S_{a}=S^{(0)}_{1}\) and \(S_{b}=(M,\dots,M)\), respectively.

We note here that the system of linear equations given by (4) has a direct correspondence with the structure of the infinitesimal generator Q above. By ordering states in levels as above, one can express (5) as

Note that elements in the diagonal of each block A_ii, for 0≤i≤N−1, are equal to zero. Moreover, the sub-vectors \(\textbf {b}_{0}^{(a)}\) and \(\textbf {b}_{1}^{(a)}\) are related to the boundary conditions in (4). In particular, each element in \(\textbf {b}_{0}^{(a)}\) (equivalently, \(\textbf {b}_{1}^{(a)}\)) is obtained by dividing the element corresponding to the same row in \(\textbf {q}^{(a)}_{0}\) (equivalently, \(\textbf {q}_{1}^{(a)}\)), which corresponds to a given state S in level \({\cal \hat S}^{(0)}\subset {\mathcal {C}}\) (equivalently, \({\mathcal {S}}^{(1)}\subset {\mathcal {C}}\)), by Δ_S.

We note that expressing (5) in the QBD structure above is relevant for solving it efficiently. In particular, instead of solving this system in terms of p^(a)=(I−A)⁻¹b^(a), one can propose a forward-elimination-backward-substitution approach, which leads to Algorithm 1.

Time until absorption

In this section, we apply similar arguments to the ones above in order to efficiently solve (7), which allows one to compute the first order moments of the time to absorption for an SVT process with one absorbing state; see the “Summary statistics” section. For convenience we assume here that the absorbing state is \(S_{a}=(M,\dots,M)\), the only state in level \({\mathcal {S}}^{(N)}\). Thus, the infinitesimal generator is given by

Given (7) and the QBD structure of the matrix Q above, and following similar arguments to those in the previous sections, one can rewrite (7) in matrix form as

In the equation above, each vector \(\textbf {m}_{j}^{(1)}\) contains first order moments \(m_{S}^{(1)}\) for states \(S\in {\mathcal {S}}^{(j)}\), for 0≤j≤N−1. According to (7), each row of each matrix \(\hat {\mathbf {A}}_{jj'}\) (for j^′∈{j−1,j,j+1}), which corresponds to a state \(S\in {\mathcal {S}}^{(j)}\phantom {\dot {i}\!}\), is obtained by dividing the same row in \(\textbf {Q}_{jj'}\phantom {\dot {i}\!}\) by Δ_S, but with elements in the diagonal equal to zero. Finally, each element of each sub-vector c_j, for 0≤j≤N−1, which corresponds to a state \(S\in {\mathcal {S}}^{(j)}\), is equal to \(\frac {1}{\Delta _{S}}\). It is clear that Algorithm 1 can be adapted for solving this system, leading to Algorithm 2.

Graph-automorphism lumping

It is easy to prove that states in the same orbit belong to the same level for an SVT process, since graph-automorphisms only permute vertex-states. That is, for any 1≤j≤r, if S_k,S_l∈L_j then there exists a 1≤p≤N such that \(S_{k},S_{l}\in {\mathcal {S}}^{(p)}\). This has a direct implication on the structure of the matrices V and W, and thus the structure of the lumped infinitesimal generator \(\hat {\mathbf {Q}}\) in (13).

Recall that the lumped infinitesimal generator \(\hat {\mathbf {Q}}\) describes the transition rates between orbits, so the particular order in which the orbits are considered has an impact on the structure of \(\hat {\mathbf {Q}}\). Let us organise orbits into groups according to the level of the states within the orbits. That is, given the partition \({\mathcal {L}}=\{L_{1},L_{2},\dots,L_{r}\}\), we consider \(\{\mathcal {L}^{(0)},\mathcal {L}^{(1)},\dots,\mathcal {L}^{(N)}\}\) where \(\mathcal {L}^{(j)}\) contains all the orbits with states in \({\mathcal {S}}^{(j)}\). For example, if there is a lumping corresponding to a partition with r=10 orbits, for N=5, where states in L₁ and L₂ belong to \({\mathcal {S}}^{(0)}\), states in L₃ belong to \({\mathcal {S}}^{(1)}\), states in L₄ and L₅ belong to \({\mathcal {S}}^{(2)}\), states in L₆, L₇ and L₈ belong to \({\mathcal {S}}^{(3)}\), states in L₉ belong to \({\mathcal {S}}^{(4)}\) and states in L₁₀ belong to \({\mathcal {S}}^{(5)}\), then

Once orbits and states are ordered as above (by levels), it is easy to check that matrices V and W are diagonal-by-blocks,

This means that the graph-automorphism lumping approach, when applied to an SVT, which is a QBD, leads to another QBD, i.e.

This also means that Algorithms 1 and 2 can be used when solving the lumped systems of equations in Section 5.1. For Algorithm 2, matrices \(\hat {\mathbf {A}}_{ij}\), for j∈{i−1,i,i+1}, are replaced by \(\textbf {W}_{i}\hat {\mathbf {A}}_{ij}\textbf {V}_{j}\), and vectors c_j are substituted by W_jc_j for all 0≤j≤N−1. In practice, it is better if the matrices \(\textbf {W}_{i}\hat {\mathbf {A}}_{ij}\textbf {V}_{j}\) can be determined without having to compute the matrix multiplication. Algorithm 1 can also be adapted, and this is particularly straight forward when the absorbing state is fixed by the automorphism group.

So far, we have presented two exact methods to analyse network dynamics. In the first case, if the model satisfies (11) then network automorphisms can be used to reduce the state-space dimension via lumping. In the second case, we have shown that the infinitesimal generator of an SVT process has the structure of a QBD, and hence efficient algorithms exist to compute summary statistics of interest. While the second case is relatively easy to check for a given model, in the first case it is generally not practical to check that all elements of the automorphism group and pairs of states satisfy (11). Thus in this section we describe some broad classes of models that satisfy (11), so that one only needs to check whether the model of interest belongs to one of these classes. This builds on recent work (Ward and Evans 2019).

Homogeneous SVT processes

Recall that in an SVT process only one vertex, the transition vertex, changes vertex-states during a single transition, thus it is only necessary to specify the transition rates between pairs of vertex-states to define an SVT process. While this rate may depend on various quantities (such as the number of neighbours of the transition vertex in each vertex-state), when the dependence of the transition rate function on these quantities is independent of the particular vertex involved, we call the model a homogeneous single-vertex transition (HomSVT) model. Consequently, the dynamics of the model can be represented graphically by a directed network on the vertex-states, where edges represent possible transitions. We call this a network dynamics diagram and an example for the SIS model of epidemic spreading is illustrated in Fig. 1. Table 1 is a list of HomSVT processes and Fig. 2 illustrates the corresponding network dynamics diagrams.

SIS model. A graphical representation of the SIS model of epidemic spreading, indicating that a susceptible (S) vertex becomes infected (I) at a rate βm_I, where m_I is the number of infected neighbours, and an infected vertex becomes susceptible at a rate γ

Full size image

HomSVT network diagrams. Network diagrams for HomSVT processes listed in Table 1. Vertex-states are indicated by labels within nodes of the networks, possible transitions by directed edges and the transition rates are written alongside the corresponding edges, rate parameters are listed in Table 1. The degree of the transition vertex is denoted by k and m_I is the number of neighbours of the transition vertex with vertex-state I, and similarly for other vertex-states, e.g. m_C, m_± etc. In (h), f(x,q)=x^q+ε[1−x^q−(1−x)^q] (see (Castellano et al. 2009) for details). The transition in (o) and (p) have rate 1 if the inequalities are satisfied

Full size image

Note that high accuracy mean-field approximations for HomSVT models have been derived (Gleeson 2011; 2013; Fennell and Gleeson 2019) and it has also been shown that the transition rates of HomSVT models can depend on any metric that is invariant under automorphisms (Ward and Evans 2019), which is natural for typical network dynamics models (otherwise the network structure would be ignored). Thus the full network automorphism group can be used to lump state-space, i.e. G=H.

We will now illustrate our analysis on a very simple numerical example of a HomSVT model where it is possible to obtain analytical results.

Biased voter model

Inspired by the voter model of opinion dynamics, we introduce a continuous-time HomSVT which we dub the “biased voter” model, and consider its dynamics on the two-edge path network illustrated in Fig. 3, which for brevity we will refer to as simply the path network.

Path network. The two-edge path network on which we consider the Biased Voter model

Full size image

Vertex states represent opinions of which there are only two, denoted by 0 or 1, and so \(\mathcal {W}=\{0,1\}\). If vertex v has degree k and opinion 0, it switches to opinion 1 at a rate αm₁/k, where α≥0 is a rate parameter and m₁ is the number of v’s neighbours that have opinion 1. Similarly, if v has opinion 1 then it switches to opinion 0 at a rate βm₀/k, where β≥0 is a rate parameter and m₀ is the number of v’s neighbours that have opinion 0. If α=β, then this model corresponds to the standard voter model in continuous time, however if α>β then there is a bias towards opinion 1. We can set α=1 and consider α≥β without loss of generality, since we can always relabel the opinions, relabel and re-scale the rates, and re-scale time.

The opinion dynamics are then described in terms of a continuous-time Markov chain \({\mathcal {X}}=\{\textbf {X}(t):\ t\geq 0\}\), where X(t) represents the opinions of the population at time t and can take values among states in the state-space \({\mathcal {S}}=\{000,001,010,100,011,101,110,111\}\), where the vertex-states of vertices 1 to 3 are indicated from left to right. The consensus states 000 and 111 are absorbing and with the ordering of states

the biased voter model dynamics on the path network are described in terms of the infinitesimal generator

However, from Fig. 3 we observe that swapping vertices 1 and 3 leaves the graph unchanged. This is the only non-trivial symmetry in the graph and thus the automorphism group is the cyclic group of order 2. Swapping vertices 1 and 3 in the state 001 results in 100, thus these states can be lumped together. Similarly, states 011 and 110 can be lumped together and all other states are fixed by this swap. Consequently, the lumping partition can be written as \({\mathcal {L}}=\{L_{1},L_{2},\dots,L_{6}\}\), where L₁={000}, L₂={001,100}, L₃={010}, L₄={011,110}, L₅={101} and L₆={111}. For this example, the collector and distributor matrices from the “Graph-automorphism lumping” section are

respectively, thus the lumped infinitesimal generator is

We now illustrate for the original and lumped CTMCs how to compute probabilities of absorption into the states 000 and 111, given by

In order to compute these probabilities, (4) leads to

Taking into account the boundary conditions \(p_{000}^{(000)}=1\) and \(p_{111}^{(000)}=0\), we get the system in matrix form

which inherits the structure of Q for the transient (i.e. non-absorbing) states. We note that in this case, because 000 is a single state within its level, b₁=A₁₀. The lumped system is

where \(\tilde {\mathbf {p}}_{i}\) is just the lumped version of p_i, and \(p_{i}^{(000)}={\lim }_{t\rightarrow \infty }\mathbb {P}(\textbf {X}(t)=000\ | \ \textbf {X}(0)\in L_{i})\). To apply Algorithm 1, we note that since 000 is the only state in the zeroth level, H₀=1 and J₀=0. Then if \(\tilde {\mathbf {A}}_{12}=\textbf {W}_{1}\textbf {A}_{12}\textbf {V}_{2}\) and \(\tilde {\mathbf {A}}_{21}=\textbf {W}_{2}\textbf {A}_{21}\textbf {V}_{1}\), we also have

Since H₂ is a two-by-two matrix, we can easily find \(\tilde {\mathbf {p}}_{2}=\textbf {H}_{2}^{-1}\textbf {J}_{2}\) and \(\tilde {\mathbf {p}}_{1}=\tilde {\mathbf {A}}_{12}\tilde {\mathbf {p}}_{2}+\textbf {J}_{1}\), resulting in

In Fig. 4, we illustrate the lumped transition rates and the probabilities of absorption for this example for the case where α=1 and β=0.5. We use a single state in each orbit to represent that orbit (e.g. the state 100 represents the lumping cell L₂), and the number of states within each orbit is indicated by the size of the grey circles. The lumped transition rates are illustrated by arrows and the size of the arrow head is scaled according to the corresponding transition rate. The coloured circles indicate the probabilities of absorption; those above each state correspond to absorption into the state 111 and those below correspond to absorption into the state 000. The colour and size of the circles corresponds to the absorption probability according to the colour bar. From Fig. 4 it is clear that states 010 and 101 are not accessible from any other states. Moreover, the bias means that it is more likely that the absorbing state 111 is reached from 010 than the absorbing state 000 is reached from 101.

Path network transitions. Lumped transition rates and probabilities of absorption for the Biased Voter model on the two-edge path network with α=1 and β=0.5. Orbits are indicated by a single representative state, e.g. 000. The size of the grey circles indicates the size of the corresponding orbit, arrows illustrate non-zero transition rates and coloured circles correspond to absorption probabilities into 000 (below the orbit representative) and 111 (above the orbit representative)

Full size image

Heterogeneous SVT processes

It is also possible to use network automorphisms to lump the state-space of heterogeneous SVT (HetSVT) models, if the heterogeneity respects (at least some of) the network symmetries. As an example, we consider an SIS epidemic spreading model inspired by a peripatetic health-care worker network (Temime et al. 2009; López-García and Kypraios 2018), illustrated in Fig. 5. This network represents the spread of a hospital-acquired infection in an intensive care unit where four patients (nodes 4, 5, 6 and 7) are treated by different health-care workers (nodes 1, 2 and 3). The different types of health care lead to different infection rates β₁, β₂ and β₃. We consider an SIS epidemic model where all nodes recover at rate γ, and where the rate at which each node becomes infected (\(S\rightarrow I\)) is the sum of all rates β_i over edges with infected neighbours. In particular, for the case illustrated in Fig. 5, if one considers the notation β₁₄=β₄₁=β₁, β₂₄=β₄₂=β₂, β₂₅=β₅₂=β₃, etc, and β_ij=0 if vertices i and j are not connected, then each susceptible vertex i becomes infected at a rate \(\sum _{j}\beta _{ji}S(j)\), where S(j)=0 if j is susceptible, and S(j)=1 if j is infective.

Peripatetic health-care worker network with heterogeneous infection rates

Full size image

The size of the state-space for this example is 2⁷=128, so we will not write out the full state-space or the infinitesimal generator. As in the Biased Voter model example, we will denote states by 0000000, 0000001 etc., where the vertex-states of vertices 1 to 7 are indicated from left to right. The state 0000000 is the only absorbing state.

We can represent this type of HetSVT model using a weighted network, where the positive weights correspond to the heterogeneous infections rates. The automorphism group of this weighted network preserves weights as well as edges and consequently satisfies (11), so we can use these symmetries to lump state-space. For the example illustrated in Fig. 5, the automorphism group is generated by a single permutation, which in cycle notation is (2 3)(4 6)(5 7), and is isomorphic to the cyclic group of order 2. Note that the automorphism group of the weighted graph is a sub-group of the automorphism group of the unweighted graph. However, we will assume the particular case where β₂=β₃. In this case, the generators of the automorphism group of the weighted network are (4 5) and (2 3)(4 6)(5 7), and this group is isomorphic to the dihedral group D₈, which is also the automorphism group of the unweighted graph. The corresponding lumped state-space consists of 42 orbits.

In Fig. 6 we illustrate the structure of the infinitesimal generator for the case γ=1.0, β₁=0.12, β₂=β₃=0.35.

Lumped transition rates for the SIS model on the peripatetic health-care worker network. Orbits are indicated by a single representative state that is formatted in correspondence with Figure 5. The size of the grey circles indicates the size of the corresponding orbit, arrows illustrate non-zero transition rates and the size of the arrow head scales with the corresponding transition rate

Full size image

We choose a single state to represent each orbit in the lumping partition and these are illustrated in correspondence with Fig. 5, so that the top vertex is vertex 1, with 2 below on the left, 3 below on the right and so on. The size of each orbit is indicated by the grey circle behind the orbit representative, and the sizes of these are either 1, 2, 4 or 8. The possible transitions between orbits are indicated by arrows between the orbit representatives. As in Fig. 4, the sizes of the arrow heads are scaled according to the corresponding transition rate. Note that for each level between 1 and 4, there is at least one state that cannot be reached by an infection, e.g. the state 0011100 in level 3.

In Fig. 7 we illustrate the mean times to absorption from each orbit in the lumping partition for the cases β₁=0.12, 1.0 and 2.0, each with γ=1.0 and β₂=β₃=0.35.

Lumped mean absorption times for the SIS model of epidemics on the peripatetic health-care worker network illustrated in Fig. 5, with recovery rate γ=1 and infection rates β₁=0.12, 1 and 2, and β₂=β₃=0.35. Initial orbits S are indicated by orbit representatives on the horizontal axis and the mean times to absorption m_S is plotted on the vertical axis

Full size image

The initial orbits are indicated on the horizontal axis by the corresponding orbit representatives and the time to absorption is plotted on the vertical axis. The orbits are grouped together by levels and these are indicated by the shading. Note that intuitively the largest absorption times are from the state 1111111, but the absorption time increases dramatically from β=1 to β=2. While the absorption times generally increase with the level, there are several cases where the absorption time is higher in a lower level state than some of the higher level states, indicating the importance of network effects.

In this paper we have proposed a unified framework to represent network dynamics in terms of continuous-time Markov chains. We have shown how a sub-set of these models (SVT processes) lead to the analysis of QBD processes in the theory of CTMCs, for which there are a number of algorithmic techniques already available in the literature. Instead of focusing on studying the master equation, our interest has been the computation of a number of summary statistics related to absorption (namely, probabilities of absorption when two absorbing states are present in the system, and time until absorption when there is just one single absorbing state). We have extended the results in (Kiss et al. 2017; Simon et al. 2011) to a broader class of models in which network symmetries can be used to reduce the size of state-space via lumping, we have also illustrated how the lumping approach can be represented in a matrix-oriented fashion, and the interplay between this and the QBD structure of SVT processes.

We believe that our results should generalise to broader classes of network, for example heterogeneous SVT processes could be re-framed as dynamics on weighted networks. Another possible future direction for this work is to consider how the lumping criterion that we have presented, (11), could be applied to more general types of model symmetry, including cases where there isn’t an underlying network. For example, there is potential to exploit lumping based on model symmetry in other disciplines including queuing theory, cellular and stochastic automata (Buchholz and Kemper 2004), Markovian agent-based models (Banisch et al. 2013; Banisch and Lima 2015; KhudaBukhsh et al. 2019), hierarchical Markov models (Buchholz 1995), continuous-time Bayesian networks (Nodelman et al. 2002) and more general CTMC models.

The interpretation of the summary statistics related to absorption depend on the process under consideration, cf. consensus in the voter model with the end of an infection in the SIS epidemic model. However, it is clear that our approach generalises to other summary statistics amenable to analysis by means of first-step arguments. This includes quasi-stationary distributions (Darroch and Seneta 1967; Gómez-Corral, A and López-García 2012), the probability distribution of the number of events (e.g., infections) during a particular time interval (e.g., until end or detection of the outbreak) (Economou et al. 2015; Gómez-Corral and López-García 2017; López-García et al. 2019), the probability of a particular node (i.e., individual) suffering a particular state (e.g., an infection) during a particular time interval (e.g., during an outbreak) (López-García 2016), the area under the curve of infectives for epidemic processes (Ball and Clancy 1995), or stationary (i.e., steady-state) distributions for non-absorbing processes (Economou et al. 2015).

Exact analyses of Markov chain dynamics on networks are rare since authors typically resort to mean-field or related approximations. These can do surprisingly well, but it is unclear what model and network features ensure this (Gleeson et al. 2012). One of the simplest mean-field approximations corresponds to the ‘well-mixed’ case, whose master equation can be viewed as the network automorphism lumping of the fully connected network (Kiss et al. 2017; Simon et al. 2011). It would be interesting to identify other graphs where symmetry gives rise to significant lumping and better dynamical approximations. For example, it has been shown that the number of lumped states is linear in the number of vertices for the complete and star graphs, and is cubic in the number of vertices for the household graph (Simon et al. 2011). Other graphs that have high symmetry, and hence a (relatively) small number of lumped states include bi- and multi-partite graphs, regular trees, graphs with multiple isomorphic disconnected components and graphs with many leaves. However, our experience is that identifying such graphs can be extremely difficult as they necessarily have to have a large amount of symmetry and this makes computing the lumped state-space difficult. An alternative approach might be to consider how lumping could be used to quantify the accuracy of mean-field approximations, for example by considering how close certain partitions of states are to exact lumping partitions, or approximate lumping based on structurally similar sets of vertices, rather than vertices that can be permuted via automorphisms, e.g. using local symmetries (KhudaBukhsh et al. 2019). There is also scope to develop more accurate mean-field approximations, similar to those in (Gleeson 2011; 2013; Fennell and Gleeson 2019), for the type of unified framework proposed here.

In summary, we have presented a general framework of dynamics on networks. One of the main purposes of this general model is to show how a wide range of dynamical processes on networks can be formulated as Markov chains. Our hope is that this will help to formalise the analytical treatment of more complex and more realistic models, ultimately improving our understanding of the phenomena that they represent.

Martín López-García acknowledges support from the Medical Research Council (UK) through the Skills Development Fellowship, MR/N014855/1, and from the Ministry of Science, Innovation and Universities (Government of Spain), through the project PGC2018-097704-B-I00

1. Jonathan A. Ward and Mart\'{i}n L\'{o}pez-Garc\'{i}a contributed equally to this work.

Authors and Affiliations

1. Department of Applied Mathematics, School of Mathematics, University of Leeds, Leeds, LS2 9JT, UK

   Jonathan A. Ward & Martín López-García

Contributions

Both authors conceived the idea, reviewed the literature, and contributed to the analysis in Sections 2-4, and to writing, reading and revising the manuscript. JAW developed the code that produced Figs. 4, 6 and 7.

Corresponding author

Correspondence to Jonathan A. Ward.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

Reprints and permissions