# Mean Field Dynamics of Graphs I: Evolution of Probabilistic Cellular
Automata for Random and Small-World Graphs

Lourens J. Waldorp<EMAIL_ADDRESS>Jolanda J. Kossakowski  Department of
Psychology
University of Amsterdam
Nieuwe Achtergracht 129-B
1018 XE Amsterdam
The Netherlands

###### Abstract

It was recently shown how graphs can be used to provide descriptions of
psychopathologies, where symptoms of, say, depression, affect each other and
certain configurations determine whether some- one could fall into a sudden
depression. To analyse changes over time and characterise possible future
behaviour is rather difficult for large graphs. We describe the dynamics of
networks using one-dimensional discrete time dynamical systems theory obtained
from a mean field approach to (elementary) probabilistic cellular automata
(PCA). Often the mean field approach is used on a regular graph (a grid or
torus) where each node has the same number of edges and the same probability
of becoming active. We show that we can use variations of the mean field of
the grid to describe the dynamics of the PCA on a random and small-world
graph. Bifurcation diagrams for the mean field of the grid, random, and small-
world graphs indicate possible phase transitions for certain parameter
settings. Extensive simulations indicate for different graph sizes (number of
nodes) that the mean field approximation is accurate. The mean field approach
allows us to provide possible explanations of ‘jumping’ behaviour in
depression.

cellular automata, discrete dynamical system, nonlinear dynamics, bifurcation

###### pacs:

64.60.aq en 05.45.-a

##  I Introduction

In psychopathology sudden changes from ’normal’ to depressed moods can occur.
Such ’discontinuities’ can be the result of a relatively small change in the
environment or person. Recently, mental disorders have been described as a
network of interacting symptoms Borsboom  _et al._ (2011); van Borkulo  _et
al._ (2014), which provides a framework where an explanation for such sudden
changes in mood could be found. A mental disorder can be viewed as a network
of symptoms, each symptom influencing other symptoms. For instance, lack of
sleep during the night could lead to poor concentration during the day, which
in turn could lead to lack of sleep again by worrying that your job may be on
the line. Here we use this idea and model the dynamics of psychopathology
networks as probabilistic cellular automata. Then to analyse the dynamics we
use a mean field approach where each node is similar in behaviour to all
others. We extend known results for the mean field approach in this context to
other types of graphs (random and small-world graphs), where the mean field
can be interpreted as a weighted average of all nodes in the graph. We also
give concentration inequalities for approximations using variations of the
grid mean field.

Cellular automata are discrete dynamical systems that have deterministic,
local rules to move from one generation to the next Wolfram (1984a); Sarkar
(2000). Introduced by Von Neumann, the most famous version is Conways game of
life, popularised by Gardner, and has found many applications from computer
science Wolfram (1984b) to neuronal population modelling Kozma  _et al._
(2005) to epdemiology Kleczkowski and Grenfell (1999). In a cellular automaton
each cell or node in a finite grid (usually a subset of $\mathbb{Z}^{2}$) can
be ’active’ or ’inactive’ (1 or 0) and if, for instance, two (direct)
neighbours are active, then the node will become active at the next time step.
Another example of a cellular automaton is bootstrap percolation, where each
node can only become active and cannot be inactivated by its neighbours, and
the objective is to determine the initial configuration of active nodes that
result in all nodes being active Janson  _et al._ (2012). In general, a new
generation in a cellular automaton is determined by a local and homogeneous
update rule $\phi$. For each node $x$ in the graph this induces a sequence of
states, an orbit. A (random) configuration at time 0 then determines whether
all nodes in the network will be active, inactive, or whether the network will
demonstrate periodic behaviour. A generalisation of a cellular automaton is to
introduce a probability $p_{\phi}$ to decide whether or not a node will become
active or not determined by a node’s neighbours. One such rule is the majority
rule which gives the probability to switch depending on whether the majority
of its neighbours are active. Such a system is called a probabilistic cellular
automaton (PCA). Here we will investigate the dynamic behaviour of the
proportion of active nodes (density) for PCA with a majority rule that are
defined on toroidal, random and small-world graphs.

Many versions of PCA exist and of particular interest are those that behave
similar to the Ising network. The reason is that the Ising network is often
used to model realistic phenomena, like magnetisation Kindermann  _et al._
(1980); Sethna  _et al._ (2004) or psychopathologies van Borkulo  _et al._
(2014). We have in mind the application to psychopathology here. In such
systems the symptoms of disorders are the nodes in the graph and edges between
the symptoms are estimated from data using the Ising mdoel van Borkulo  _et
al._ (2014) or from verbal accounts. In our companion paper in this issue we
elaborate on real data analysis using results from discrete Markov chains.

Watts showed that a one-dimensional, large-scale cellular automaton
(deterministic) where the connectivity between nodes was arranged as a small-
world, could perform the density (all zeroes) and synchronisation (alternating
all zeroes and all ones) tasks. Newman and Watts gave approximations for path
length and clustering on a small-world, to obtain an analytic solution to the
threshold above which a large number of nodes are active. Callaway  _et al._
also studied percolation in different graph topologies in deter, focussing on
the consequences of (randomly) deleting nodes. Here again the objective was to
concentrate on stable solutions of the graphs. In a probabilistic version,
Tomassini  _et al._ investigated a one-dimensional PCA on a regular and small-
world graph in terms of its performance on the density and synchronisation
tasks. They determined by using evolutionary algorithms that a small-world
topology is most efficient to solving both tasks, corresponding to the results
of Watts in a deterministic version. Their objective was different from ours
in that here we are interested in all types of dynamic behaviour (stable or
not), and specifically representing this behaviour for the PCA by the mean
field.

Our starting point is the work by Balister  _et al._ and Kozma  _et al._ where
a two-dimensional (toroidal) grid on a PCA is defined. The mean field is then
used to determine the unconditional probability distribution of the density
(relative number of active nodes). Balister  _et al._ show that the mean field
model predicts a bifurcation for small values of the probability of a node
switching to another state and determine its critical point for a
neighbourhood of size five (see also Kozma  _et al._ , 2004, 2005). This is of
particular interest in our case as it may explain mood disorders (e.g.,
depression or manic-depression) from symptoms and their connectivity. To apply
these results to random and small-world graphs we determine the marginal
distribution across the possible degree probabilities given the topology of
the random or small-world graph. Extending results of homogeneous graphs has
been applied to social networks Barrat  _et al._ (2008) and to cellular
automata Janson  _et al._ (2015).

We first introduce probabilistic cellular automata in Section II. Then in
Section III we show how the traditional version of a PCA on a grid can be
reduced to a single discrete time dynamical system, called the mean field. In
Section III.2 we show that for the random graph we can use a variation on the
formulation for the grid of the dynamical system to describe dynamics. We use
these results on the random graph to show in Section III.3 that we can obtain
a similar approximation for the small-world graph, again using the formulation
for the grid. Having shown that these approximations are appropriate, we see
in Section IV what the dynamics of the process is for the different
topologies. We follow these theoretical results by extensive simulations to
verify the accuracy of the mean field in Section V.

##  II Probabilistic cellular automata

A cellular automaton is a dynamical system of nodes in a fixed, finite grid
where directly connected nodes determine the state of a node at each
subsequent time step Wolfram (1984b). Each node $x$ in a node set
$V=\\{1,2,\ldots,n\\}$ is at time $t$ in one of the states of a finite
alphabet $\Sigma$. The nearest neighbours in the graph $G=(V,E)$ are given by
the edges in $E$. Often the graph $G$ is the square lattice $\mathbb{Z}^{2}$,
where each node has exactly four neighbours (see e.g., Grimmett, 2010). A
local rule determines based on the direct neighbours what the value of the
alphabet of node $x\in V$ will be at time $t+1$. Let the neighbourhood of $x$
be the set of nodes that are directly connected to $x$, $\Gamma(x)=\\{y\in
V:(x,y)\in E\\}$. A local rule $\phi:\Gamma\to\Sigma$ assigns for each
configuration of the neighbourhood of $x$ a value $a\in\Sigma$. If we
additionally introduce a probability for each local configuration $\phi$, we
obtain a probabilistic cellular automaton (PCA). The probability is a function
$p_{a}:\Sigma^{|\Gamma|}\times\Sigma\to[0,1]$ such that a probability is
assigned to each node $x$ to have label $a$ for a configuration $\phi$
dependent on the neighbourhood $\Gamma(x)$, with $\sum_{a\in\Sigma}p_{a}=1$.
The local rule $\phi$ is applied iteratively to each result, and hence induces
a stochastic process with sequence $\Phi_{t}:V\to\Sigma$ for each time step
$t$. For node $x\in V$ we write $\Phi_{t}(x)=\phi^{t}(x_{0})$, where $x_{0}$
is the value at time $t=0$, and $\Phi_{t}(V)$ is the image for all nodes
simultaneously. Each node therefore has an orbit Hirsch  _et al._ (2004),
which is the sequence $(\Phi_{t}(x),t\geq
0)=(x_{0},\phi(x_{0}),\phi(\phi(x_{0})),\ldots)$.

  * Example 1

(Majority rule on 0/1) Let the alphabet be $\Sigma=\\{0,1\\}$ and take a
finite subset of the square lattice $V\subset\mathbb{Z}^{2}$. In this lattice
$V$ each node has 4 (nearest) neighbours. We use the majority rule, which says
that if $|\phi_{x}^{-1}(1)|=|\\{y\in\Gamma(x):\Phi_{t}(y)=1\\}|$ is greater
than $2$, then the node $x$ will be 1 with probability $1-p$, and otherwise 0
with probability $p$. The sequence $(\Phi_{t}(x),t\geq 0)$ is any orbit of 0s
and 1s $(0,1,1,\ldots)$.

In Example II the probability $p$ determined by the neighbourhood $\Gamma(x)$
is independent of the state of the node $x$ itself. Such a model is called
totalisitc Balister  _et al._ (2006). Additionally, the probabilities for 0
and 1 were defined by the same parameter $p$, which is then called symmetric,
i.e., $p_{1}=1-p_{0}$. Here we focus on totalistic and symmetric models with
size 2 alphabet $\Sigma=\\{0,1\\}$.

##  III Mean field approximation on graphs

The key ingredient of the mean field approximation, shown by Balister  _et
al._ , is that the properties of interest are uniform over the graph. For the
(toroidal) grid topology this is easy to see: Any node $x\in V$ has the same
number of neighbours $|\Gamma|=4$, where each node in the neighbourhood
becomes 0 or 1 by the same local rule. It follows that any four nodes in the
grid could serve as part of the neighbourhood for $x$. In a probabilistic
automaton, therefore, the local rule depends only on the number of 1s in any
random draw of $4$ nodes from all nodes $V$. We first consider the case for a
grid and then move onto the random and small-world graph.

###  III.1 Mean field on a grid

Let the graph $G_{\rm grid}(n,\Gamma)$ be a grid with $n$ nodes and boundary
conditions such that each node has exactly four neighbours. We consider the
density $\rho_{t}$ defined by $|\Phi_{t}^{-1}(1)|/n$, where the set of nodes
that are 1 is $\Phi_{t}^{-1}(1)=\\{y\in V:\phi^{t}(y)=1\\}$. It follows that
we require the probability of $\\{\Phi_{t}(x)=1\\}$ given a certain number of
nodes in state 1 at the previous time point. The probability of switching to
state 1, $\xi_{|\Gamma|}(r)$, is conditional on $r$ of the neighbours that are
$1$. Let $|\phi^{-1}(1)|=|\\{y\in\Gamma:\phi(y)=1\\}|$ be the number of 1s in
the neighbourhood of $x$. Then we define the probability of state 1 given that
$r$ neighbours are 1 as

$\displaystyle\xi_{|\Gamma|}(r)=\mathbb{P}(\Phi_{t}(x)=1\mid|\phi^{-1}(1)|=r)$
(1)

One of the possibilities to define $\xi_{|\Gamma|}$ is the majority rule: if
the majority of neighbours in $\Gamma$ are 1, then the node will be 1 with
probability $1-p$ at $t+1$. The majority rule is defined as

$\displaystyle\xi_{|\Gamma|}(r)=\begin{cases}p&\text{if }r\leq|\Gamma|/2\\\
1-p&\text{if }r>|\Gamma|/2\end{cases}$ (2)

To obtain the probability of $\Phi_{t}(x)=1$, the state of $x$ being in state
1, we need to determine the probability of a neighbourhood having $r$ active
nodes. In the mean field we have that the probability of a 1 is homogeneous
and so we obtain a binomial distribution for the number of 1s in the
neighbourhood $\Gamma$. The probability that the number of nodes in the
neighbourhood equals $r$ is in the mean field the same as $|\Gamma|$ Bernoulli
trials each with probability of success $\rho_{t}$. Hence

$\displaystyle\mathbb{P}(|\phi^{-1}(1)|=r\mid\rho_{t})=\binom{|\Gamma|}{r}\rho_{t}^{r}(1-\rho_{t})^{|\Gamma|-r}$
(3)

where $r=0,1,\ldots,|\Gamma|$. Then the probability of the event
$\\{\Phi_{t}(x)=1\mid\rho_{t}\\}$, that node $x$ will have state 1 at time
$t+1$ given the density $\rho_{t}$ at time $t$, is

$\displaystyle p_{\rm
grid}(\rho_{t})=\sum_{r=0}^{|\Gamma|}\xi_{|\Gamma|}(r)\binom{|\Gamma|}{r}\rho_{t}^{r}(1-\rho_{t})^{|\Gamma|-r}.$
(4)

Let

$\displaystyle
q_{|\Gamma|/2}(\rho_{t})=\sum_{r=0}^{|\Gamma|/2}\binom{|\Gamma|}{r}\rho_{t}^{r}(1-\rho_{t})^{|\Gamma|-r}$

Combining this probability with the majority rule (2) gives

$\displaystyle p_{\rm
grid}(\rho_{t})=pq_{|\Gamma|/2}(\rho_{t})+(1-p)\left(1-q_{|\Gamma|/2}(\rho_{t})\right)$

which is used in Kozma  _et al._. This mean field result follows directly from
Theorem 2.1 in Balister  _et al._. Let $B(n,p)$ denote a binomial random
variable with $n$ Bernoulli trials each with success probability $p$.

  * Lemma 2

(Balister  _et al._ , 2006, Theorem 2.1) Let $G_{\rm grid}(n,\Gamma)$ be a
grid with a PCA as defined above with $\xi_{|\Gamma|}(r)$ according to the
majority rule in (2). Then the evolution of the number of active nodes
$n\rho_{t}$ is

$\displaystyle n\rho_{t+1}=B(n,p_{\rm grid}(\rho_{t})).$ (5)

The mean and variance for the density $\rho_{t}$, respectively, $\mu_{\rm
grid}=p_{\rm grid}$ and $\sigma^{2}_{\rm grid}=p_{\rm grid}(1-p_{\rm
grid})/n$.

It follows that the conditional probability of $n\rho_{t+1}$ active nodes in
the grid on $\\{n\rho_{t}=k\\}$ at the previous time step is

$\displaystyle\mathbb{P}(n\rho_{t+1}=r\mid n\rho_{t}=k)=$
$\displaystyle\qquad\binom{n}{r}p_{\rm grid}(\rho_{t})^{r}(1-p_{\rm
grid}(\rho_{t}))^{n-r}$ (6)

It is easily seen that this is a discrete time Markov process on a finite
state space of size $n$, since the probability of $n\rho_{t+1}$ depends only
on $\rho_{t}$. Equation (III.1) is the transition probability of the discrete
time Markov process for the number of active nodes in the graph.

Because $\phi_{t}(x)$ is Bernoulli distributed $B(1,p_{\rm grid}(\rho_{t}))$
for all nodes $x\in V$, and the number of active nodes $|\Phi^{-1}(1)|$ is the
sum of these Bernoulli trials, we can apply the law of large numbers so that
for large $n$, $\rho_{t+1}$ is close to $\mu_{t}:=\mu(\rho_{t})$ with high
probability. Indeed, we can use Chernov’s bound to suggest that using $p_{\rm
grid}$ is good enough for large graphs.

  * Lemma 3

(Accuracy bound of density) Let $n\rho_{t}=\sum_{x\in V}\phi^{t}(x)$ be the
sum of $n$ Bernoulli trials given by (5), with mean of the density $p_{\rm
grid}(\rho_{t})$. For every $0<\varepsilon<\min\\{p_{{\rm grid}},1-p_{\rm
grid}\\}$, let $\delta=2\exp(-\varepsilon^{2}/2\sigma^{2}_{{\rm grid}})$. We
then have with probability at least $1-\delta$

$\displaystyle|\rho-p_{{\rm grid}}|\leq\sqrt{\frac{p_{\rm grid}(1-p_{\rm
grid})}{n}2\log(2/\delta)}$ (7)

A proof is in the Appendix. So we can use the mean field $p_{\rm
grid}(\rho_{t})$ for grids of large size $n$. With $\delta=0.05$, we obtain
the interval with probability at least 0.95 of $[\mu_{\rm
grid}-2.72\sigma_{\rm grid},\mu_{\rm grid}-2.72\sigma_{\rm grid}]$. Another
interval can be obtained from the DeMoivre-Laplace central limit theorem. This
theorem tells us that for large enough $n$, $z_{\rm grid}=(\rho_{t}-p_{{\rm
grid}})/\sigma_{{\rm grid}}$ is distributed as $N(0,1)$. In fact, if the third
order moment of $z_{\rm grid}$ is $c<\infty$, the Berry-Esseen theorem says
that the order of approximation of the distribution of $\rho_{t}$ to the
normal distribution is $O(3c/\sqrt{n})$ Venkatesh (2013). This provides an
interval for $\rho_{t+1}$ as a measure of accuracy with $[\mu_{{\rm
grid}}-1.96\sigma_{{\rm grid}},\mu_{{\rm grid}}+1.96\sigma_{{\rm grid}}]$ with
probability $0.95$. Clearly, in both limit laws the size of the network $n$
determines the accuracy of the approximation.

###  III.2 Mean field on a random graph

In the original setting of a grid (with boundary conditions, so a torus) the
number of neighbours is fixed and it was seen that the mean field
approximation $p_{\rm grid}$ was accurate for the density because each node is
identical with respect to a change depending on its neighbours. Here we
introduce the neighbourhood size $|\Gamma|$ as a random variable and then
determine the probability of $\Phi_{t}(x)=1$ given $\rho_{t}$ by averaging
over all possible sizes of neighbourhoods weighted by its probability for
neighbourhood size. This is done in a random graph where each node has a
binomial number of neighbours. Let $G_{\rm rg}(n,p_{e})$ be a random graph
with $n$ nodes and (constant) probability $p_{e}$ of an edge being present
Bollobás (2001); Durrett (2007). Let the size of the neighbourhood $|\Gamma|$
be a binomial random variable with maximal value $n-1$ neighbours and
probability $p_{e}$, that is, $B(n-1,p_{e})$ and

$\displaystyle\mathbb{P}_{p_{e}}(|\Gamma|=k)=\binom{n-1}{k}p_{e}^{k}(1-p_{e})^{n-k-1}.$

Then the probability of obtaining an active node can be defined conditionally
on the event $\\{|\Gamma|=k\\}$, the neighbourhood having size $k$. Then
marginalising over the possible the neighbourhood size, we obtain $p_{\rm rg}$
for the probability of a node being active in the binomial process. Proofs can
be found in the Appendix.

  * Lemma 4

(Probability on a random graph) Consider a PCA on a random graph $G_{\rm
rg}(n,p_{e})$ with edge probability $p_{e}$ and a local rule $\phi$ determined
by the majority rule for the neighbourhood $|\Gamma|=k$ according to $\xi_{k}$
as in (4). Then the probability of obtaining an active node at time $t+1$ is

$\displaystyle p_{\rm rg}(\rho_{t})=$
$\displaystyle\sum_{k=0}^{n-1}\sum_{r=0}^{k}\xi_{k}(r)\binom{k}{r}\rho_{t}^{r}(1-\rho_{t})^{k-r}\mathbb{P}_{p_{e}}(|\Gamma|=k).$
(8)

Intuitively we would expect that if we restrict the neighbourhoods in (Lemma
4) to the expected neighbourhood size $p_{e}(n-1)$ for each node, then the
approximation should be reasonably close. This would make it possible to
analytically determine fixed points more easily and simplify computation for
large graphs considerably. We next show that this is a reasonable approach.

  * Lemma 5

(Simple probability on a random graph) If the neighbourhood in $G_{\rm
rg}(n,p_{e})$ of each node in a PCA as in Lemma III.2 is fixed with the
expected number of nodes under the random graph $\nu=\lfloor
p_{e}(n-1)\rfloor$ such that $k=\nu$ in $p_{\rm grid}$, then the probability
$p_{\rm rg}$ reduces to

$p_{\rm
grid}^{\nu}(\rho_{t})=\sum_{r=0}^{\nu}\xi_{\nu}(r)\binom{\nu}{r}\rho_{t}^{r}(1-\rho_{t})^{\nu-r}.$
(9)

The approximation error is

$\displaystyle|p_{\rm rg}-p_{\rm grid}^{\nu}|\leq|p-1/2|$ $\displaystyle\times
2\exp(-(n-1)\varepsilon^{2}/p_{e}(1-p_{e})+\log(n)),$

as $\varepsilon$ decreases to 0 and with $0<p_{e}<1$.

This implies of course that the number of active nodes $n\rho_{t}$ in the
random graph with probability $p_{\rm rg}$ in (Lemma 4) and the number of
active nodes with probability $p_{\rm grid}^{\nu}$ in (9) converge in
probability with exponential rate with graph size $n$.

Figure 1: The expectation $p_{\rm rg}$ (blue, solid curve) of equation (Lemma
4) and $p_{rmrand}$ (red, dotted curve) of equation (9) with $p=0.1$ and
$p_{e}=0.3$. Left panel shows the curves for a graph of size $n=25$, showing a
clear difference between the curves, and the right panel for graph size
$n=100$. Note that the difference between the curves at the crossings with the
$45^{\circ}$ line is small.

  * Remark

To retain a probability of an edge in $p_{\rm grid}^{\nu}$ leads to a larger
approximation error, i.e., using

$\displaystyle p_{\rm grid}^{n-1}(\rho_{t})=$
$\displaystyle\sum_{r=0}^{n-1}\xi_{\nu}(r)\binom{n-1}{r}(\rho_{t}p_{e})^{r}(1-\rho_{t}p_{e})^{n-r-1}$
(10)

leads to an error of at most $|p-1/2|$, which makes it non ignorable (see the
Appendix).

We now have expression (9) similar to (4) for a random graph with the
probability of an active node at time $t$ determined by both the density
$\rho_{t}$ and an edge being present $p_{e}$ in the size of the neighbourhood.
From (Lemma 4) in Lemma III.2 and Corollary III.1 and Lemma III.2 the
evolution equation for the random graph follows.

  * Theorem 7

(Evolution on a random graph) Let a PCA with local rule $\phi$ be defined on a
random graph $G_{\rm rg}(n,p_{e})$, and let the probability of obtaining an
active node in $G_{\rm rg}(n,p_{e})$ be $p_{\rm grid}^{\nu}$ as defined in
(9). Then the evolution equation for the random graph with large $n$ is

$\displaystyle n\rho_{t+1}=B(n,p_{\rm grid}^{\nu}(\rho_{t})).$ (11)

The mean and variance of the density $\rho_{t}$ respectively, $\mu_{\rm
grid}^{\nu}=p_{\rm grid}^{\nu}$ and $\sigma_{\rm grid,\nu}^{2}=p_{\rm
grid}^{\nu}(1-p_{\rm grid}^{\nu})/n$.

We immediately have that the probability $p_{\rm grid}^{\nu}$ is close to the
density $\rho$ for each time point $t$ for large graph size $n$. In fact, we
find by the triangle inequality

$|\rho-p_{\rm grid}^{\nu}|\leq|\rho-p_{\rm rg}|+|p_{\rm rg}-p_{\rm
grid}^{\nu}|,$

and both terms converge to 0. The first term $|\rho-p_{\rm rg}|$ converges to
0 by Lemma III.1 with $p_{\rm rg}$ and using Lemma III.2, and $|p_{\rm
rg}-p_{\rm grid}^{\nu}|$ converges to 0 by Lemma III.2.

The process $n\rho_{t}$ on a random graph is also a discrete time Markov
process, as before, and has transition probability (III.1) with $p_{\rm
grid}^{\nu}$. We can then apply a similar analysis of dynamics to $\mu_{\rm
grid}^{\nu}=p_{\rm grid}^{\nu}$ as before.

Note that we require for obvious reasons that the graph is connected. It
follows that we need a minimum probability $p_{e}$ such that the graph is
connected. The probability that a random graph $G_{\rm rg}$ is connected is
$\exp(-\exp(-\lambda))$, where $p_{e}=(\log n+\lambda+o(1))/n$ with $\lambda$
fixed (Bollobás, 2001, Theorem 7.3). For instance, if we choose the
probability of $G_{\rm rg}$ being connected to be 0.99 and we use $n=50$, then
we obtain $\lambda=4.6$ and hence $p_{e}=0.17$. We can therefore not go below
0.17 for a graph with $n=50$ nodes.

###  III.3 Small-world graph

A small-world graph is one which has high average clustering and low average
path length, relative to a random graph with the same number of nodes and
edges. These graphs have been shown to model realistic networks like those of
working relations between actors and the nerve cells in the worm C. elegans
(Watts and Strogatz, 1998), and subsequently the small-world has been shown to
apply to many different networks, like the (parcellated) brain (Sporns and
Honey, 2006). And most recently, the network of symptoms as defined by the
diagnostic statistical manual (a compendium to diagnose patients) has been
found to be a small-world. This finding is a possible explanation for the
correlations between pairs of symptoms found in different subpopulations
(Borsboom  _et al._ , 2011).

Here we use the modified Newman-Watts (NW) small-world of Newman and Watts
(1999), where for a given grid structure where each node has neighbourhood
$\Gamma$, a set of $(n-1)p_{w}$ edges is on average independently added to the
graph, where $p_{w}$ is the probability of two nodes being wired. Such a graph
is denoted by $G_{\rm sw}(n,\Gamma,p_{w})$. The same idea as with the random
graph, where the probability for an active node was corrected by the
probability of the degree of a node, averaged over all possible neighbourhood
sizes, can be used for the random part in the NW small-world. In the NW small-
world we start with a grid with neighbourhood size $|\Gamma|$, which is fixed,
and augment the graph randomly with edges according to a binomial variable
with probability $p_{w}$. Let

$\displaystyle\mathbb{P}_{p_{w}}(|\Gamma|=k)=\binom{n-1}{k}p_{w}^{k}(1-p_{w})^{n-k-1}$

We then obtain

$\displaystyle p_{\rm sw}(\rho_{t})=$
$\displaystyle\sum_{k=|\Gamma|}^{n-1}\sum_{r=0}^{k-|\Gamma|}\xi_{k}(r)\binom{k}{r}\rho_{t}^{r}(1-\rho_{t})^{k-r}\mathbb{P}_{p_{w}}(|\Gamma|=k).$
(12)

We could define the small-world probability using this definition. But we can
split up $p_{\rm sw}$ in two terms, one involving the fixed neighbourhood
$\Gamma$ of the grid, and one random neighbourhood consisting of the possible
shortcuts. We therefore start with the probability in a grid $p_{\rm grid}$
corrected by the $(1-p_{w})^{n-|\Gamma|}$ requiring that no possible randomly
added edges are present, i.e., we obtain

$p_{\rm grid}^{\rm sw}=p_{\rm grid}(1-p_{w})^{n-|\Gamma|}$ (13)

for the first part of the fixed grid. Then, in accordance with the random part
of the NW small-world, a probability is added to emulate the possible
additional neighbours in the random part of the graph, ignoring the first
$|\Gamma|$ neighbours from the grid. Define the probability

$\displaystyle p_{\rm rg,\Gamma}(\rho_{t})=$
$\displaystyle\sum_{k=|\Gamma|+1}^{n-1}\sum_{r=0}^{k-|\Gamma|}\xi_{k}(r)\binom{k}{r}\rho_{t}^{r}(1-\rho_{t})^{k-r}\mathbb{P}_{p_{w}}(|\Gamma|=k),$
(14)

where the first $|\Gamma|$ neighbours are ignored since they were included
already as neighbours in the grid structure in $p_{\rm grid}^{\rm sw}$. Then
we can write the small-world probability as $p_{\rm sw}=p_{\rm
grid}(1-p_{w})^{n-|\Gamma|}+p_{\rm rg,\Gamma}$. For the second part, however,
we have the approximation as before from the random graph, leaving out the
first $\Gamma$ nodes from the grid. This leads to the simplification using the
grid probability only

$p_{{\rm
grid},\Gamma}^{\nu}(\rho_{t})=\sum_{r=|\Gamma|+1}^{\nu}\xi_{\nu}(r)\binom{\nu}{r}(\rho_{t})^{r}(1-\rho_{t})^{\nu-r},$
(15)

where $\nu=\lfloor p_{w}(n-|\Gamma|)\rfloor$. The error of approximation using
$p_{{\rm grid},\Gamma}^{\nu}$ instead of $p_{{\rm rg},\Gamma}$ follows
immediately from Lemma III.2 for fixed grid neighbourhood $\Gamma$, except
that the first $|\Gamma|$ nodes in the grid are taken out.

  * Corollary 8

Let $G_{\rm sw}(n,\Gamma,p_{w})$ be the NW small-world graph of size $n$ with
$|\Gamma|$ nodes in the fixed neighbourhood for each node. Furthermore, let
$0<p_{w}<1$ be the wiring probability and $\nu=\lfloor
p_{w}(n-|\Gamma|)\rfloor$. Then the approximation error for the probability
using the grid structure $p_{\rm grid}^{\nu}$ in (15) in the random part is

$\displaystyle|p_{\rm rg,\Gamma}-p_{\rm grid,\Gamma}^{\nu}|\leq|p-1/2|$
$\displaystyle\times 2\exp(-(n-|\Gamma|)\varepsilon^{2}/p_{w}(1-p_{w})$
$\displaystyle\quad+\log(n-|\Gamma|+1)),$

for $\varepsilon>0$.

Equations (III.3) to (15) and Corollary III.3 prove the following equation for
the evolution on an NW small-world.

  * Theorem 9

(Evolution on a small-world) Define a PCA on a Newman-Watts small-world graph
$G_{\rm sw}$ of size $n$ with $|\Gamma|$ neighbours for the initial graph and
wiring probability $p_{w}$. Then with probabilities $p_{\rm grid}^{\rm sw}$ in
(13) and $p_{\rm grid,\Gamma}^{\nu}$ in (15), the evolution of the number of
active nodes is

$\displaystyle n\rho_{t+1}$ $\displaystyle=B(|\Gamma|,p_{\rm grid}^{\rm
sw}(\rho_{t}))$ $\displaystyle\quad+B(n-|\Gamma|,p_{{\rm
grid},\Gamma}^{\nu}(\rho_{t}))$ (16)

with mean and variance for the density $\rho_{t}$, respectively,

$\mu_{\rm sw}=\frac{|\Gamma|}{n}p_{\rm grid}^{\rm
sw}+\frac{n-|\Gamma|}{n}p_{{\rm grid},\Gamma}^{\nu}$ (17)

and

$\displaystyle\sigma^{2}_{\rm sw}$ $\displaystyle=\frac{|\Gamma|}{n^{2}}p_{\rm
grid}^{\rm sw}(1-p_{\rm grid}^{\rm sw})$
$\displaystyle\quad+\frac{n-|\Gamma|}{n^{2}}p_{{\rm
grid},\Gamma}^{\nu}(1-p_{{\rm grid},\Gamma}^{\nu}).$ (18)

We write $p_{\rm sw}^{\nu}=p_{\rm grid}^{\rm sw}(|\Gamma|/n)+p_{\rm
grid,\Gamma}^{\nu}((n-|\Gamma|)/n)$ for the NW small-world probability based
on the approximation with the grid. Figure 2 shows two examples of the
approximation $p_{\rm grid,\Gamma}^{\nu}$ for the random part in the NW small-
world. It is clear from the corollary that convergence is a bit slow for small
graphs since the difference of nodes in the fixed neighbourhood $\Gamma$ and
in the expected $p_{w}(n-|\Gamma|)$ neighbours in the random part, determines
the rate. Again, we can use $p_{\rm sw}^{\nu}$ to determine the dynamics of
the mean field.

Figure 2: The expectation $p_{\rm sw}$ (blue, solid curve) of equation
(LABEL:eq:pr-sw-rg) and $p_{\rm sw}^{\nu}$ (red, dotted curve) of equation (9)
with $p=0.1$ and $p_{w}=0.3$. Left panel shows the curves for a graph of size
$n=25$, showing a clear difference between the curves, and the right panel for
graph size $n=100$. Note that the difference between the curves at the
crossings with the $45^{\circ}$ line is small.

##  IV Dynamics of the mean field

To investigate the dynamics we treat the mean field function $p_{\rm grid}$
for the grid, $p_{\rm grid}^{\nu}$ for the random graph, and $p_{\rm sw}$ for
the NW small-world as a discrete dynamical system. We can then determine in
principle the fixed points and describe its behaviour in the long term. In
general, however, obtaining the fixed points is not trivial. Indeed, Balister
_et al._ (2006) provide an analytical solution for the fixed points for a
specific case in $G_{\rm grid}$, but mention that a general solution is
difficult. Janson  _et al._ (2015) give analytical solutions for the fixed
points when leaving out the majority rule, making it a deterministic system.
Here we keep the majority rule sacrificing the possibility of determining the
critical points analytically. We therefore describe the qualitative behaviour
of $p_{\rm grid}$, $p_{\rm grid}^{\nu}$, and $p_{\rm sw}$.

###  IV.1 Dynamics of the mean field in a grid

The dynamics of the mean field in the grid $G_{\rm grid}$ from (4) have been
described in Balister  _et al._ (2006) and Kozma  _et al._ (2005) for a
neighbourhood size of $|\Gamma|=5$. The function $\mu_{\rm grid}=p_{\rm grid}$
is continuous and since $[0,1]$ is closed and bounded, we find that $\mu_{\rm
grid}$ has at least one fixed point in $[0,1]$ (Holmgren, 1996; Hirsch  _et
al._ , 2004). A fixed point is one where we find $p_{\rm
grid}(\rho_{t})=\rho_{t}$. Finding the fixed points for $p_{\rm grid}$ in
general is not trivial. Balister  _et al._ (2006) showed that if $|\Gamma|=5$
in the finite grid, then $p=7/30\approx 0.233$ is a critical point, such that
if $p$ is in $[7/30,1/2]$ then there is a stable fixed point at $\rho=0.5$,
but when $p<7/30$ then $\rho=0.5$ is unstable and there are two other stable
fixed points. This can be seen in Figure 3, which shows two bifurcation plots,
where for each value of $0<p\leq 0.5$ the function $\mu_{\rm grid}=p_{\rm
grid}$ is iteratively applied for about 1000 steps, and only the last 50 are
plotted. Figure 3 shows that for $|\Gamma|=5$ neighbours the fixed point at is
stable for $p\in[0,7/30)$ and bistable for $p\in[7/30,0.5]$, as predicted.

Figure 3: Bifurcation plots of $p_{\rm grid}$ in a graph of size $n=100$ for
$|\Gamma|=5$ and $|\Gamma|=15$ neighbours.

Since $p_{\rm grid}$ is continuous, stability can be checked by considering
the derivative $\partial\mu_{\rm grid}/\partial\rho=\dot{\mu}_{\rm grid}$. If
$|\dot{\mu}_{\rm grid}|$ is bounded by 1, then the fixed point $\rho$ is
attractive, otherwise it is repellent. The derivative with respect to
$\rho_{t}$ is

$\displaystyle\dot{\mu}_{\rm grid}(\rho_{t})=$
$\displaystyle\sum_{r=0}^{|\Gamma|}\xi_{|\Gamma|}(r)\binom{|\Gamma|}{r}(r-\rho_{t}|\Gamma|)\rho_{t}^{r-1}(1-\rho_{t})^{|\Gamma|-r-1}$

For example, the derivative for $p=0.15$ is not bounded by 1 for all values of
$\rho$; the fixed point $\rho=0.5$ is repellent since at this point
$\dot{\mu}(0.5)\approx 1.359$, and so iteration of $\mu$ will lead away from
0.5. The derivative for $p=0.35$ is smaller than 1 (0.672) and so in this case
$\rho=0.5$ is an attractive fixed point. It can be seen that for $|\Gamma|=5$
the critical point is at 0.233, as predicted by theory (Balister  _et al._ ,
2006). It can also be seen that increasing the neighbourhood size to
$|\Gamma|=15$ (right panel) increases the critical point. This increase in
critical point corresponds to the simulations in Kozma  _et al._ (2005) where
(’long range’) edges were added to the nodes, which increased the
neighbourhood size.

###  IV.2 Dynamics of the mean field in a random graph

The dynamics of $p_{\rm grid}^{\nu}$ in the random graph $G_{\rm rg}$ are
similar to that of the grid. The main difference is that the critical point of
the bifurcation is closer to $p=1/2$. As is clear from the definition of
$p_{\rm grid}^{\nu}$ in (9), the only difference with that of the grid is the
neighbourhood size which is increased to $\nu=\lfloor p_{e}(n-1)\rfloor$.
Figure 4 shows the result for a graph with $n=25$ nodes (left panel) and for a
graph with $n=100$ nodes. The approximation of $p_{\rm grid}^{\nu}$ is quite
accurate, also for the location of the critical point. With a graph of size
$n=100$ the accuracy is such that $p_{\rm rg}$ and $p_{\rm grid}^{\nu}$ are
nearly indistinguishable, which corresponds to the result in Lemma III.2.

Figure 4: Plots of $\mu_{\rm rg}$ (red) and $\mu_{\rm grid}^{\nu}$ (blue) as a
function of $p$. In the left panel the bifurcation plots are given for a
random graph of size $n=25$ and in the right panel for a random graph of size
$n=100$. All plots are obtained with edge probability $p_{e}=0.4$.

###  IV.3 Dynamics of the mean field in a small-world

The dynamic behaviour of $p_{\rm sw}$ is shown in Figure 5. Generally, the
behaviour is similar to that on the random graph. In Figure 5 the left panel
shows a bifurcation plot of $p_{\rm sw}$ and and $p_{\rm sw}^{\nu}$ on $G_{\rm
sw}(49,0.4)$. The accuracy of $p_{\rm sw}^{\nu}$ improves greatly for larger
$n$, as seen in the right panel of Figure 5 for $G_{\rm sw}(100,0.4)$.

Figure 5: Bifurcation plots of the small-world mean field $\mu_{\rm
sw}^{\nu}(\rho)$ based on the grid with $\nu=\lfloor p_{e}(n-1)\rfloor$
neighbours (blue) and the mean field $\mu_{\rm sw}(\rho)$ based on the all
possible neighbours in the random graph (red). In the left panel a small-world
of $n=25$ nodes and in the right panel a graph of $n=100$ nodes; all graphs
are obtained with the probability of wiring (adding edges) in the NW small-
world of $p_{w}=0.4$.

For low values of new edges in the NW small-world $p_{w}$, the probability
$p_{\rm sw}^{\nu}$ is smaller than in the grid. This is because the
probability $p_{\rm grid}^{\rm sw}=p_{\rm grid}(1-p_{w})^{n-|\Gamma|}$ is
corrected by the number of edges not added to the graph.

##  V Numerical evaluation of the mean field

To evaluate the accuracy of the predictions of the mean field in the grid,
random, and small-world graph, we simulated networks of different sizes in the
topology of a grid, a random graph, and a small-world graph. For each
simulation run 0/1 data were generated according to one of the three types of
graph using the R package IsingSampler van Borkulo  _et al._ (2014). The in
combination with the majority rule the PCA was run for a certain duration $T$
and the average states of the last section of the time series was determined
to see if it matches that of the predictions of the mean field. To determine
the accuracy we used both $90\%$ and $95\%$ confidence intervals obtained from
the central limit theorem for each of the three different graphs (see Section
III.1).

For each combination of parameters, 100 graphs were simulated in the topology
of an unweighted grid, a random graph, and a small-world graph. We varied the
size of the network $n\in\\{16,25,49,100\\}$, the number of time points
$T\in\\{50,100,200,500,5000\\}$, and the probability of an active node in the
majority rule $p\in\\{0.1,0.2,0.3,0.4,0.5\\}$, see (2). We also varied the
probability of an edge in the random graph $p_{e}\in\\{0.1,0.2,\ldots,0.9\\}$,
and the probability of wiring in the small-world graph
$p_{w}\in\\{0.1,0.2,\ldots,0.9\\}$. For $t=0$, a random number of nodes was
set to active by using the R package _IsingSampler_ version 0.2 Epskamp
(2013). Figures 6 and 7 show the topology of some simulated random graphs (6)
and small-world graphs (7) arranged in a square.

$n=16$$p_{e}=0.1$$n=25$$n=49$$n=100$

$p_{e}=0.2$

$p_{e}=0.3$

$p_{e}=0.4$

$p_{e}=0.5$

Figure 6: Examples of simulated random graphs with network size
$n\in\left\\{16,25,49,100\right\\}$ and
$p_{e}\in\left\\{0.1,0.2,0.3,0.4,0.5\right\\}$.

All simulated data, figures, as well as the used R-code are publicly available
at the OSF Asante  _et al._ (2016).

$n=16$$p_{w}=0.1$$n=25$$n=49$$n=100$

$p_{w}=0.2$

$p_{w}=0.3$

$p_{w}=0.4$

$p_{w}=0.5$

Figure 7: Examples of simulated small world graphs with network size
$n\in\left\\{16,25,49,100\right\\}$ and
$p_{w}\in\left\\{0.1,0.2,0.3,0.4,0.5\right\\}$.

Results in bifurcation diagrams with 90% and 95% confidence intervals are
shown in Figure 8. Mean density estimates, visualized with red dots, were
calculated by dividing each simulation into snippets in which the distance
between density estimates ($\delta$) did not exceed $0.4$. Other values for
$\delta$ were considered, but it was found that changing $\delta$ to $0.3$ or
$0.5$ resulted in similar percentages. We calculated for each snippet the
percentage of mean density estimates that fell within a 90% and 95% confidence
interval. Table 1 shows the mean percentages, marginalized over the number of
time points $T$, the variation in $p_{e}$, $p_{w}$ and the number of nodes $n$
obtained from the central limit theorem with bounds $\hat{\rho}\pm 1.96$ for
the 95% confidence interval, and $\hat{\rho}\pm 1.64$ for the 90% confidence
interval. Figure 9 gives a 3-dimensional representation of the percentage of
mean density estimates that fall within a 95% confidence interval. Results
from the 90% confidence interval are not presented, as these were similar to
the results from the 95% confidence interval. It can be seen that the mean
field approximation accurately estimates the density of the network structures
across various simulation conditions.

As seen in Figure 9, a local dip occurs for all network structures at $p=0.3$
and becomes more extreme as $n$ increases in size. Also shown in Figure 8, the
mean density estimates at $p=0.3$ fall less often in the 95% confidence
interval in comparison to the mean density estimates at other values for $p$.
This is partly due to the fact that the standard error, a parameter needed for
the calculation of the confidence intervals, depends on the network size; as
the network size $n$ increases, the standard error becomes smaller as well as
the resulting confidence interval. Furthermore, as this phenomenon occurs in
all simulated network structures, we believe that this behaviour results from
the fact that the mean field approximation has a bit of trouble adjusting to
the one-phase stability, after being in an area where phase transitions may
occur. All in all, results show that the mean field approximation also
performs well when non-regular network structures are under consideration.

$t=50$$n=16$Torus$t=500$$n=49$$t=5000$$n=100$

Randomgraph

Small-worldgraph

Figure 8: Bifurcation diagrams of a torus (upper panel), a random graph
(middle panel; $p_{e}=0.5$) and a small world graph (lower panel;
$p_{w}=0.5$). Grey solid area = 90% confidence interval around bifurcation.
Dashed grey lines = 95% confidence interval around bifurcation. Red dots =
mean density estimates at different values of $p$.

$n=16$Torus$n=25$$n=49$$n=100$

Randomgraph

Small-worldgraph

Figure 9: Visualization of the percentage of mean density estimates that fall within a 95% confidence interval in 2D (torus; left column) and 3D (random graph; middle column, and small world graph; right column) at $t=500$. Table 1: Percentages of density estimates with a 90% or a 95% confidence interval. $p_{e}$ = edge probability (random graph). $p_{w}$ = rewiring probability (small-world graph). Structure |  T |  N |  $p_{e}$/$p_{w}$ |  $p$ |  90%CI |  95%CI
---|---|---|---|---|---|---
Torus |  50 |  16 |  - |  0.1 |  0.75 |  0.80
|  |  |  |  0.2 |  0.99 |  0.99
|  |  |  |  0.3 |  0.92 |  0.99
|  |  |  |  0.4 |  0.94 |  0.99
|  |  |  |  0.5 |  0.96 |  0.99
|  500 |  49 |  - |  0.1 |  0.29 |  0.67
|  |  |  |  0.2 |  0.97 |  0.99
|  |  |  |  0.3 |  0.63 |  0.75
|  |  |  |  0.4 |  0.75 |  0.83
|  |  |  |  0.5 |  0.81 |  0.87
|  5000 |  100 |  - |  0.1 |  0.12 |  0.55
|  |  |  |  0.2 |  1.00 |  1.00
|  |  |  |  0.3 |  0.50 |  0.50
|  |  |  |  0.4 |  0.75 |  0.83
|  |  |  |  0.5 |  0.75 |  0.75
Random graph |  50 |  16 |  0.5 |  0.1 |  0.97 |  0.97
|  |  |  |  0.2 |  1.00 |  1.00
|  |  |  |  0.3 |  0.55 |  0.76
|  |  |  |  0.4 |  0.78 |  0.88
|  |  |  |  0.5 |  0.85 |  0.92
|  500 |  49 |  0.5 |  0.1 |  1.00 |  1.00
|  |  |  |  0.2 |  1.00 |  1.00
|  |  |  |  0.3 |  0.98 |  0.98
|  |  |  |  0.4 |  0.92 |  0.97
|  |  |  |  0.5 |  0.94 |  0.98
|  5000 |  100 |  0.5 |  0.1 |  1.00 |  1.00
|  |  |  |  0.2 |  1.00 |  1.00
|  |  |  |  0.3 |  0.87 |  0.91
|  |  |  |  0.4 |  0.98 |  1.00
|  |  |  |  0.5 |  0.99 |  1.00
Small-world graph |  50 |  16 |  0.5 |  0.1 |  0.34 |  0.76
|  |  |  |  0.2 |  0.20 |  0.43
|  |  |  |  0.3 |  0.38 |  0.59
|  |  |  |  0.4 |  0.54 |  0.69
|  |  |  |  0.5 |  0.63 |  0.75
|  500 |  49 |  0.5 |  0.1 |  1.00 |  1.00
|  |  |  |  0.2 |  1.00 |  1.00
|  |  |  |  0.3 |  1.00 |  1.00
|  |  |  |  0.4 |  0.98 |  0.99
|  |  |  |  0.5 |  0.99 |  1.00
|  5000 |  100 |  0.5 |  0.1 |  1.00 |  1.00
|  |  |  |  0.2 |  1.00 |  1.00
|  |  |  |  0.3 |  0.99 |  0.99
|  |  |  |  0.4 |  1.00 |  1.00
|  |  |  |  0.5 |  1.00 |  1.00
Figure 10: Examples of the evolution of a torus (upper panel; $p=0.1$), a
random graph (middle panel; $p=0.3,p_{e}=0.6$) and a small world graph (lower
panel; $p=0.1,p_{w}=0.9$).

Figure 10 visualizes the evolution of selected simulation conditions. Phase
transitions were obtained for the random graph, specifically at
$p_{e}=\left\\{0.3,0.4\right\\}$ and $p=\left\\{0.3,0.4\right\\}$. We also
obtained phase transitions for the small world graph at
$p_{w}=\left\\{0.9\right\\}$ and $p=\left\\{0.1\right\\}$. For the duration in
our simulations, we did not obtain phase transitions in any other simulations.
However, as we observed stable chains throughout the simulation study, we
expect to obtain phase transitions in all conditions when $t$ approaches
infinity.

##  VI Conclusions and discussion

To model the complex dynamics of large-scale networks (graphs) is in general
difficult. This is because there are many different ’agents’ that operate
within the graph. In particular, if the nodes in the graph represent symptoms
and the edges represent their mutual influence, then the interacting symptoms
show complex behaviour on a macroscopic scale, e.g., at the level of the
number of active symptoms. Here we showed that the mean field model for a
probabilistic cellular automaton with majority rule, can serve as an accurate
approximation to such large-scale graphs, and can simplify analysis of the
dynamics of such systems. Specifically, we showed that averaging across the
different possible degrees for a random and small-world graph, results in
approximations that lie with high probability close to a generalised version
of the mean field on a torus. These theoretical results were confirmed by
extensive simulations, showing correspondence between the mean field and
simulated graph activity for different size graphs.

Our approximation is based on the formulation of the grid (torus) where a
relatively simple sum over possible active nodes determines the probability of
a randomly selected node in the graph being active. We showed that for large
graphs this approximation is accurate. This simplification could serve to
obtain a more extensive analysis of the dynamics such as that presented in
Janson  _et al._ (2015). There the majority rule (the probabilistic element)
was removed from the model, to obtain exact fixed points for the model. Here
we chose not to remove the probabilistic element since we aim to introduce
different rules for updates than the majority rule, like a conditional Ising
probability.

Our initial motivation for these results was to obtain a model where we could
assess the risk of a single person based on the estimate of the graph and the
corresponding probability of an active node $p$, to determine the risk of that
person ’jumping’ from one state into another. This risk assessment might be
useful in a clinical setting where a decision in a particular type of
intervention is required. This idea is pursued in the companion paper in this
issue.

## Appendix

###### Proof.

(Lemma III.1) Let the Kullback-Leibler divergence between $p+\varepsilon$ and
$p$ be defined as a function of $0\leq\varepsilon\leq p$

$\displaystyle h_{+}(\varepsilon)$
$\displaystyle=(p+\varepsilon)\log\frac{p+\varepsilon}{p}$
$\displaystyle+(1-p-\varepsilon)\log\frac{1-p-\varepsilon}{1-p}$

and similarly, define $h_{-}(\varepsilon)=h_{+}(-\varepsilon)$. Then Chernov’s
bound (Lesigne, 2005; Venkatesh, 2013) for the density $\rho_{t}$ of a grid
with $n$ nodes and its mean at time $t$, $p_{{\rm grid}}(\rho_{t})$ defined in
(4), for $0<\varepsilon<\min\\{p_{\rm grid},1-p_{\rm grid}\\}$ immediately
gives

$\displaystyle\mathbb{P}(|\rho_{t}-p_{{\rm grid}}(\rho_{t})|>\varepsilon)\leq$
$\displaystyle\exp(-nh_{+}(\varepsilon))+\exp(-nh_{-}(\varepsilon))$

The Kullback-Leibler divergence can be approximated quadratically by

$\displaystyle
h_{+}(\varepsilon)=\frac{\varepsilon^{2}}{2p(1-p)}+O(\varepsilon^{3})\quad\text{as
}\varepsilon\to 0$

Similarly for $h_{-}(\varepsilon)$ gives

$\displaystyle\mathbb{P}(|\rho_{t}-p_{{\rm grid}}(\rho_{t})|>\varepsilon)\leq
2\exp(-\varepsilon^{2}/2\sigma^{2}_{{\rm grid}}(\rho_{t}))$ (19)

where $\sigma^{2}_{\rm grid}=p_{\rm grid}(1-p_{\rm grid})/n$. Let
$\delta=2\exp(-\varepsilon^{2}/2\sigma^{2}_{{\rm grid}}(\rho_{t}))$ such that
$\varepsilon=\sqrt{2\sigma^{2}_{\rm grid}\log(2/\delta)}$. Then we obtain the
result with probability at least $1-\delta$. ∎

###### Proof.

(Lemma III.2) By Theorem 2.1 in Balister  _et al._ (2006) we have that the
probability of state 1 at time $t+1$ on the event $\\{|\Gamma|=k\\}$ and given
$\rho_{t}$ is

$\displaystyle\mathbb{P}(\Phi_{t}(x)=1\mid|\Gamma|=k,\rho_{t})=$
$\displaystyle\sum_{r=0}^{k}\xi_{k}(r)\binom{k}{r}\rho_{t}^{r}(1-\rho_{t})^{k-r}$

and in the random graph of Erdös-Renyi the probability of $k$ neighbours for
any node is

$\displaystyle\mathbb{P}_{p_{e}}(|\Gamma|=k)=\binom{n-1}{k}p_{e}^{k}(1-p_{e})^{n-k-1}$

It follows that the marginal
$\sum_{k}\mathbb{P}(\Phi_{t}(x)=1\mid|\Gamma|=k,\rho_{t})\mathbb{P}(|\Gamma|=k)$
is

$\displaystyle p_{\rm rg}(\rho_{t})=$
$\displaystyle\sum_{k=0}^{n-1}\sum_{r=0}^{k}\xi_{k}(r)\binom{k}{r}\rho_{t}^{r}(1-\rho_{t})^{k-r}\mathbb{P}_{p_{e}}(|\Gamma|=k)$

as claimed. ∎

###### Proof.

(Equation Remark) To obtain (Remark) assume a fixed value $\nu$ for all $k$,
$\xi_{k}=\xi_{\nu}$ for all $k$. Then $\xi_{\nu}(r)$ only depends on $r$.
First note that

$\displaystyle\binom{k}{r}\binom{n-1}{k}=\binom{n-1}{r}\binom{n-r-1}{k-r}$

Second, observe that the order of the two sums in $p_{\Phi,G}$ can be switched
since the first sum is for the sequence $(k=0,1,\ldots,n-1)$ and the second is
$(r=0,1,\ldots,k)$, which can be seen as the lower triangular of the two-
dimensional array for $(k=0,1,\ldots,n-1)$ and $(r=0,1,\ldots,n-1)$. Upon
switching to $(r=0,1,\ldots,n-1)$ and $(k=r,r+1,\ldots,n-1)$, we obtain the
upper triangle of the array. Second, by changing the order of summation and
reordering the sums, we get

$\displaystyle\sum_{r=0}^{n-1}\xi_{\nu}(r)\binom{n-1}{r}\rho_{t}^{r}p_{e}^{r}\sum_{k=r}^{n-r-1}\binom{n-r-1}{k-r}$
$\displaystyle\times(p_{e}(1-\rho_{t}))^{k-r}(1-p_{e})^{n-k-1}$

In the sum on the right we can use the binomial theorem with $m=k-r$ and
$N=n-r-1$, which gives

$\displaystyle\sum_{m=0}^{N}\binom{N}{m}(p_{e}(1-\rho_{t}))^{m}(1-p_{e})^{N-m}=$
$\displaystyle\quad(p_{e}(1-\rho_{t})+1-p_{e})^{N}$

which leads to (Remark).

For the approximation error, write $\xi_{k}(r)=p\mathbbm{1}\\{r\leq
k/2\\}+(1-p)\mathbbm{1}\\{r>k/2\\}$ and recall that $\nu$ is fixed. Then

$\displaystyle p_{\rm rg}(\rho_{t})-p_{\rm rand}(\rho_{t})=$
$\displaystyle\quad\sum_{k=0}^{n-1}\sum_{r=0}^{k}\binom{k}{r}\rho_{t}^{r}(1-\rho_{t})^{k-r}\binom{n-1}{k}$
$\displaystyle\quad\times p_{e}^{k}(1-p_{e})^{n-k-1}[\xi_{k}(r)-\xi_{\nu}(r)]$

Using Hölder’s inequality with the $\ell_{\infty}$ and $\ell_{1}$ norms, gives

$\displaystyle|p_{\rm rg}(\rho_{t})-p_{\rm rand}(\rho_{t})|$
$\displaystyle\leq$
$\displaystyle\sum_{k=0}^{n-1}\sum_{r=0}^{k}\binom{k}{r}\rho_{t}^{r}(1-\rho_{t})^{k-r}\binom{n-1}{k}$
$\displaystyle\times
p_{e}^{k}(1-p_{e})^{n-k-1}\max_{r,k}|\xi_{k}(r)-\xi_{\nu}(r)|$

The binomial theorem for the first term of the right hand side gives

$\displaystyle\sum_{k=0}^{n-1}\sum_{r=0}^{k}\binom{k}{r}\rho_{t}^{r}(1-\rho_{t})^{k-r}\binom{n-1}{k}p_{e}^{k}(1-p_{e})^{n-k-1}$
$\displaystyle=$
$\displaystyle\sum_{r=0}^{n-1}\binom{n-1}{r}(\rho_{t}p_{e})^{r}(1-\rho_{t}p_{e})^{n-r-1}=1.$

For each $r,k$ such that $r\leq k$ we have that

$\displaystyle\xi_{k}(r)-\xi_{\nu}(r)=$ $\displaystyle\quad
p(\mathbbm{1}\\{r\leq k/2\\}-\mathbbm{1}\\{r\leq\nu/2\\})+$
$\displaystyle\quad(1-p)(\mathbbm{1}\\{r>k/2\\}-\mathbbm{1}\\{r>\nu/2\\})$

The term $|\xi_{k}(r)-\xi_{\nu}(r)|$ is at most $2p-1$ if $\nu<k$ or $1-2p$ if
$\nu\geq k$ for any $r,k$, which gives the size of the error bound. ∎

###### Proof.

(Lemma III.2) If we fix $\nu=\lfloor p_{e}(n-1)\rfloor$, the expectation of
the random variable for each node of the possible number of neighbours
$B(n-1,p_{e})$, such that each $k=\nu$ in the part for the density we obtain

$\displaystyle\sum_{k=0}^{n-1}\binom{n-1}{k}p_{e}^{k}(1-p_{e})^{n-k-1}$
$\displaystyle\times\left(\sum_{r=0}^{\nu}\binom{\nu}{r}\xi_{\nu}(r)\rho_{t}^{r}(1-\rho_{t})^{\nu-r}\right),$

from which we obtain $p_{\rm grid}^{\nu}(\rho_{t})$. The approximation error
for the probabilities is then

$\displaystyle|p_{\rm rg}(\rho_{t})-p_{\rm rand}^{\nu}(\rho_{t})|=$
$\displaystyle\left|\sum_{k=0}^{n-1}(p_{\rm grid}^{k}(\rho_{t})-p_{\rm
grid}^{\nu}(\rho_{t}))\mathbb{P}_{p_{e}}(|\Gamma|=k)\right|.$

The probability of obtaining a neighbourhood size $k$ close to the expected
number of neighbours $\nu=p_{e}(n-1)$ can be obtained from the Chernov bound
in Lemma III.1, giving $\mathbb{P}(|k-\nu|\leq t)\geq
1-2\exp(-(n-1)\varepsilon^{2}/p_{e}(1-p_{e}))$, for $\varepsilon\searrow 0$.
This leads to the difference being bound by

$\displaystyle|p_{\rm rg}(\rho_{t})-p_{\rm rand}^{\nu}(\rho_{t})|\leq$
$\displaystyle\biggl{|}\sum_{k=0}^{n-1}(p_{\rm grid}^{k}(\rho_{t})-p_{\rm
grid}^{\nu}(\rho_{t}))$ $\displaystyle\times
2\exp(-(n-1)\varepsilon^{2}/p_{e}(1-p_{e}))\biggr{|}.$

Using Hölder’s inequality with the sup and $\ell_{1}$ norms, we find that the
above is

$\displaystyle\leq$ $\displaystyle\max_{k}|p_{\rm grid}^{k}(\rho_{t})-p_{\rm
grid}^{\nu}(\rho_{t})|$
$\displaystyle\times\sum_{k=0}^{n-1}2\exp(-(n-1)\varepsilon^{2}/p_{e}(1-p_{e})).$

The difference $p_{\rm grid}^{k}(\rho_{t})-p_{\rm grid}^{\nu}(\rho_{t})$ is
determined by the mismatch between $k$ and $\nu$ and is at most $2p-1$ if
$\nu<k$ and $1-2p$ if $\nu\geq k$ for any $r,k$. And so we obtain

$\displaystyle|p_{\rm rg}(\rho_{t})-p_{\rm grid}^{\nu}(\rho_{t})|\leq|p-1/2|$
$\displaystyle\times 2\exp(-(n-1)\varepsilon^{2}/p_{e}(1-p_{e})+\log(n)),$

completing the proof. ∎

###### Proof.

(Theorem III.2) In the random graph $G_{\rm rg}(n,p_{e})$ each node has the
same number $n-1$ of possible neighbours. Hence, for each node the local rule
$\phi$ is determined by $r$ out of $n-1$ possible neighbours being both
connected and 1. Lemma III.2 shows that the probability of this binomial
process is then $p_{\rm rg}$ in (Lemma 4). Then by Lemma III.2 for large $n$
the probability $p_{\rm grid}^{\nu}$ converges to $p_{\rm rg}$, and by
consequence, the number of active nodes in both processes converges. ∎

## References

  * Borsboom  _et al._ (2011) D. Borsboom, A. O. J. Cramer, V. D. Schmittmann, S. Epskamp, and L. J. Waldorp, PLoS ONE 6, e27407 (2011).
  * van Borkulo  _et al._ (2014) C. D. van Borkulo, D. Borsboom, S. Epskamp, T. F. Blanken, L. Boschloo, R. A. Schoevers, and L. J. Waldorp, Scientific reports 4 (2014).
  * Wolfram (1984a) S. Wolfram, Nature 311, 419 (1984a).
  * Sarkar (2000) P. Sarkar, ACM Computing Surveys (CSUR) 32, 80 (2000).
  * Von Neumann (1951) J. Von Neumann, Cerebral mechanisms in behavior , 1 (1951).
  * Gardner (1970) M. Gardner, Scientific American 223, 120 (1970).
  * Wolfram (1984b) S. Wolfram, Communications in Mathematical Physics 96, 15 (1984b).
  * Kozma  _et al._ (2005) R. Kozma, M. Puljic, P. Balister, B. Bollobás, and W. J. Freeman, Biological Cybernetics 92, 367 (2005).
  * Kleczkowski and Grenfell (1999) A. Kleczkowski and B. T. Grenfell, Physica A: Statistical Mechanics and its Applications 274, 355 (1999).
  * Janson  _et al._ (2012) S. Janson, T. Łuczak, T. Turova, T. Vallier, _et al._ , The Annals of Applied Probability 22, 1989 (2012).
  * Kindermann  _et al._ (1980) R. Kindermann, J. L. Snell, _et al._ ,  _Markov random fields and their applications_ , Vol. 1 (American Mathematical Society Providence, RI, 1980).
  * Sethna  _et al._ (2004) J. P. Sethna, K. A. Dahmen, and O. Perkovic, arXiv preprint cond-mat/0406320 (2004).
  * Watts (1999) D. J. Watts,  _Small-worlds: From order to chaos_ (Princeton University Press, 1999).
  * Newman and Watts (1999) M. E. Newman and D. J. Watts, Physical Review E 60, 7332 (1999).
  * Callaway  _et al._ (2000) D. S. Callaway, M. E. Newman, S. H. Strogatz, and D. J. Watts, Physical review letters 85, 5468 (2000).
  * Tomassini  _et al._ (2005) M. Tomassini, M. Giacobini, and C. Darabos, Complex Systems 15, 261 (2005).
  * Balister  _et al._ (2006) P. Balister, B. Bollobás, and R. Kozma, Random Structures & Algorithms 29, 399 (2006).
  * Kozma  _et al._ (2004) R. Kozma, M. Puljic, P. Balister, B. Bollobas, and W. Freeman, in  _Cellular Automata_, Lecture Notes in Computer Science, Vol. 3305, edited by P. M. Sloot, B. Chopard, and A. G. Hoekstra (Springer Berlin Heidelberg, 2004) pp. 435–443.
  * Barrat  _et al._ (2008) A. Barrat, M. Barthelemy, and A. Vespignani,  _Dynamical processes on complex networks_ (Cambridge University Press, 2008).
  * Janson  _et al._ (2015) S. Janson, R. Kozma, M. Ruszinkó, and Y. Sokolov, arXiv preprint arXiv:1507.07997 (2015).
  * Grimmett (2010) G. Grimmett,  _Probability on graphs: Random processes on graphs and lattices_ (Cambridge University Press, 2010).
  * Hirsch  _et al._ (2004) M. Hirsch, S. Smale, and R. Devaney,  _Differential equations, dynamical systems, and an introduction to chaos_ , 2nd ed. (Academic Press, 2004).
  * Venkatesh (2013) S. Venkatesh,  _The theory of probability_ (Cambridge University Press, 2013).
  * Bollobás (2001) B. Bollobás,  _Random Graphs_ (Springer, 2001).
  * Durrett (2007) R. Durrett,  _Random graph dynamics_ (Cambridge University Press, 2007).
  * Watts and Strogatz (1998) D. J. Watts and S. H. Strogatz, Nature 393, 440 (1998).
  * Sporns and Honey (2006) O. Sporns and C. J. Honey, Proceedings of the National Academy of Sciences USA 103, 19219 (2006).
  * Holmgren (1996) R. Holmgren,  _A first course in discrete dynamical systems_ (Springer Science & Business Media, 1996).
  * Epskamp (2013) S. Epskamp, “Isingsampler: Sampling methods and distribution functions for the ising model,” url: github.com/SachaEpskamp/IsingSampler (2013).
  * Asante  _et al._ (2016) K. Asante, E. Barbour, L. Barker, M. Benjamin, S. D. Bowman, A. P. Boughton, and L. … Chen, “Open science framework,” url: osf.io/4znzp (2016).
  * Lesigne (2005) E. Lesigne,  _Heads or tails: An introduction to limit theorems in probability_ (American Mathematical Society, 2005).