## Abstract

Applications such as neuroscience, telecommunication, online social networking, transport and retail trading give rise to connectivity patterns that change over time. In this work, we address the resulting need for network models and computational algorithms that deal with dynamic links. We introduce a new class of evolving range-dependent random graphs that gives a tractable framework for modelling and simulation. We develop a spectral algorithm for calibrating a set of edge ranges from a sequence of network snapshots and give a proof of principle illustration on some neuroscience data. We also show how the model can be used computationally and analytically to investigate the scenario where an evolutionary process, such as an epidemic, takes place on an evolving network. This allows us to study the cumulative effect of two distinct types of dynamics.

## 1. Introduction

The last decade has seen a huge rise in interest in complex networks and their applications to mass communication, social and natural phenomena. Until recently, one characteristic of such networks that is fundamental within almost all applications has received only marginal attention: that the networks may evolve (Saramäki & Kaski 2005; Kao *et al.* 2007; Borgnat *et al.* 2008; Gautreau *et al.* 2009; Vernon & Keeling 2009). This issue is distinct from network growth, or aggregative phenomena: it embodies the property that all edges within the network are transient to some extent. Very recently, general classes of dynamic networks have been proposed and studied in the theoretical computer science literature (Avin *et al.* 2008; Clementi *et al.*2008, 2009) from a complexity theory perspective. Our work looks at complementary issues driven by the need for practical tools in modelling, calibration and data analysis.

In Borgnat *et al.* (2008), it is pointed out that network evolution can be approached from distinct directions: at the macro-level, it may be observed within real datasets by studying the time course of global parameters from one snapshot to another, or at the micro-level, dynamic properties could be ascribed to the individual birth–death rules for each edge; and that specific applications will require some analytical methods to approach an *inverse problem*: given some data from a time-dependent evolving network, how best may one represent it within a suitably defined class of models?

In this paper, we consider these problems and offer an operational approach to each.

In §2, we consider micro-to-macro-models that produce evolving networks, considered as classes of Markov processes defined over the set of all possible undirected graphs on a finite set of *n* vertices. This space grows as *O*(2^{n2}), so we suggest simplifying to the case of independent dynamics for each edge.

As a further simplification, we may then allow the dynamics to depend only on the *range* of the edge. This extends the static concepts of ‘lattices plus shortcuts’ (Kleinberg 2000*a*; Newman *et al.* 2000; Higham 2007) and more general range-dependent random graphs (Grindrod 2002; Higham 2005) to the dynamic setting. In many applications, vertices are fixed within some Euclidean or underlying metric space and each edge has a natural ‘range’ representing the distance between the end vertices, which may impact on the possibility of that edge arising. For example, in traditional acquaintanceship networks physical neighbours have a good chance of knowing each other. These ideas are pervasive in the literature (e.g. the recent treatment in Franceschetti & Meester 2008). But in many applications, there is no obvious, fixed lattice topology. For example:

in communication networks, individuals may be mobile,

in online social networks, or telecommunications networks, physical position does not determine the cost or likelihood of interaction,

in cyberspace, hyperlinks are not constrained by the underlying Internet connectivity structure,

in high-frequency functional connectivity networks from neuroscience, cognitive processing tasks may be distributed across the brain, and

in proteomic interaction networks, connections may be caused by a combination of features, including electrostatic, hydrophobic or chemical similarities rather than any obvious sequence-level or geographical commonalities.

In all these cases, the concept of ‘range’ is more elusive, not corresponding directly to a simple distance, but can still have some meaning. The range of an edge reflects the transitive nature of the connection. If vertex *a* is connected to vertex *b* which is in turn connected to vertex *c*, how likely is it that vertex *a* is also connected to vertex *c*? If it is very likely then we will say that connection is (locally) transitive and the edge from *a* to *c* is short range, if it is very unlikely then the edge from *a* to *c* is long range. Cliques are full of short-range edges. Hence, when we are presented with data from an evolving graph as just a transient set of connections between an arbitrarily ordered list of vertices, there is potential to add insight by inferring a range for every possible edge. This is the *inverse problem* we will address: that of representing a given evolving graph as a range-dependent evolving graph within an imposed class and thus inferring a range for every possible edge.

Section 2 briefly discusses evolving networks in a general setting. In §3, we review the formal concept of range dependency, introduce the new extension to evolving networks and look at the inverse problem of discovering the ranges. The spectral algorithm that we develop is then illustrated in §4 and tested on some real brain activity networks. In §§5 and 6, we turn our attention to propagation problems defined on evolving graphs. Unlike percolation problems for static graphs, or dynamic propagation of infections or information on static graphs, here the dynamics and the sequential behaviour of the evolving graph interferes with the dynamics assumed for propagation of changes to vertex-dependent states. We thus consider the role of evolving graphs in transmission and threshold behaviour as well as achievable expected path lengths and shortest transit times. We believe that this work taps into an exciting and novel field for analysis with many potential applications—from information theory and neuroscience through to viral, or buzz, marketing.

## 2. Independent edge-dependent dynamics

Let *V*_{n} denote a set of *n* labelled vertices. Let *S*_{n} denote the space of all undirected graphs defined on *V*_{n}. Then |*S*_{n}|=2^{n(n−1)/2}. Any element of *S*_{n}, denoted by *a*, may be represented by an *n*×*n* symmetric adjacency matrix, which we will denote by *A*.

We will consider discrete-time Markov processes defined over *S*_{n}, whose paths consist of time-dependent sequences of elements, {*a*_{j}} in *S*_{n}, each representing the *evolving graph* at time *t*_{j}=*j* *δ**t*. Even when we restrict to the case where the transition matrix is time independent, to fully specify such a process, in general, requires 2^{n(n−1)} non-negative graph-to-graph transition probabilities. We therefore continue with a simplified class of models where the time-dependent appearance or disappearance of each individual edge is governed by a random process that is independent of all other edges. More precisely, consider an evolving graph {*a*_{j}} defined as follows. Let *α*(*e*) denote the probability that any edge, *e*, not part of the network at time *t* may be added to it over the time step *δ**t*. Let *ω*(*e*) denote the probability that any edge, *e*, that is part of the network at time *t* will be removed over the time step *δ**t*. So *α* and *ω* specify the birth and death probabilities, respectively, that we assume to be *O*(*δ**t*) for *δ**t* small.

For any pair *a*, *a*′∈*S*_{n}, let *P*(*a*′|*a*) denote the probability that *a*_{j+1}=*a*′ given that *a*_{j}=*a*, and let *E*(*a*) denote the set of edges belonging to *a*. Then, the ‘independent edge-dependent’ model yields the graph-to-graph transition probability
2.1

This expression gives the probability that exactly the right subsets of edges are added and deleted to achieve the required transition. As a result of the edge independence assumption, there are now *n*(*n*−1) parameters (*α*(*e*) and *ω*(*e*)) rather than the 2^{n(n−1)} required in the general case.

It is straightforward to calibrate such a model, given sufficient data. Suppose that we observe a sequence . Then we may estimate the parameters *α*(*e*) and *ω*(*e*) independently for each edge, using Laplace’s law from the sequence data. Specifically, suppose that an edge *e* is absent within the first *J*−1 terms of the observed sequence of graphs on exactly *M*(*e*) occasions and of these graphs exactly *m*(*e*) are followed by graphs that contain *e* (so *e* appears in those transitions), then we have the following estimate
2.2
Similarly, suppose an edge *e* occurs within the first *J*−1 terms of the observed sequence of graphs on exactly *M*^{★}(*e*) occasions and of these graphs, exactly *m*^{★}(*e*) are followed by graphs that do not contain *e* (so *e* is lost in the transition). Then
2.3

We may then use these estimates in equation (2.1) to generate any graphto-graph transition probability that is required. In particular, we could simulate and analyse sequences of networks.

## 3. Range dependency

### (a) Review of the static case

In Grindrod (2002), the class of *range-dependent random graphs* was introduced as a parameterized model that can reproduce important properties seen in real networks. Protein–protein interaction data were used to motivate and justify the concept. Closely related models based on similar principles include:

the original small world networks of Watts & Strogatz (1998), and their counterparts based on adding shortcuts rather than rewiring (Newman

*et al.*2000), where edges are either long or short range,the two-dimensional lattice-based model of Kleinberg (2000

*a*,*b*), andthe geometric model used by Pržulj and co-workers (Pržulj

*et al.*2004; Kuchaiev*et al.*2009) to describe protein–protein interactions.

Range-dependent random graphs are best introduced by imagining the vertices set out in a line and labelled by their integer positions. To simplify things further, if *n*, the number of vertices, approaches infinity, then we may approximate the graph by one on infinitely many vertices (…,−2,−1,0,1,2,…), since the edge effects become less important.

Range-dependent random graphs are then defined as follows: an edge is present between any vertices *i*_{1} and *i*_{2} with probability *p*_{i1,i2}=*f*(|*i*_{1}−*i*_{2}|), where *f* is a given monotonically decreasing function of the edge *range* |*i*_{1}−*i*_{2}|. Thus, in this model the presence or absence of an edge depends only on its range, and each edge is independent. Letting *P*_{k} be the consequent probability that any vertex has degree *k*, the generating function can be used to study the Watts–Strogatz clustering coefficient (Watts & Strogatz 1998) and the mean degree (see Grindrod 2002 for details).

In the protein–protein interaction case, and in most realistic network scenarios, the vertices will be labelled in a way that does not reflect any edge range information. Therefore, there is a natural inverse problem of reordering the vertices to reveal the range dependency. Genetic algorithms (Grindrod 2002) and more efficient spectral methods (Higham 2003; Grindrod *et al.* submitted) have proved successful in this context.

Our aim is now to develop these static ideas into the evolving graph framework.

### (b) Evolving range-dependent random graphs

Consider a set of *n* vertices labelled by location *i*=1,…,*n*. As in §2 we consider a discrete Markov process over *S*_{n} where all edges evolve independently. Suppose further that each edge *e* has transition probabilities that depend only on its range: if an edge connects vertices *i*_{1} and *i*_{2}, then we will write *k*(*e*)=|*i*_{1}−*i*_{2}| to denote its range. Then, an evolving range-dependent random graph has birth and death transition probabilities
given by functions *f*_{α}(*k*) and *f*_{ω}(*k*) that map the positive integers onto [0,1].

Now let *p*(*e*,*j*) denote the probability that the edge *e* is present within *a*_{j}, the graph at time *t*_{j}. Then using the transition probabilities above, we have the dynamical equation

A steady distribution must then satisfy
3.1
which depends only upon *k*(*e*). Hence, at equilibrium any single observation of the evolving graph appears as a range-dependent random graph with each edge present according to this range-dependent probability function *p*_{0}(*k*).

Now consider the natural inverse problem. Given an observed sequence, , with the vertices in some given ordering *i*=1,…,*n*, how can we best represent that evolving graph within the class of evolving range-dependent random graphs? It is clearly reasonable to reorder the vertices with a mapping *q*(*i*) so as to maximize the likelihood of the actual observations.

The simplest way forward is to consider each edge in turn within the evolving sequence. Suppose that *e*=(*i*_{1},*i*_{2}), in the original ordering, is observed on exactly *r*_{i1,i2} occasions (and is absent on *J*−*r*_{i1,i2} occasions). Let *R* denote the symmetric non-negative matrix with elements *r*_{i1,i2}. Under a reordering *q*(*i*) the edge range becomes *k*(*e*)=|*q*(*i*_{1})−*q*(*i*_{2})| and the likelihood of the observations for this edge is given by
Trivially, we can rewrite this as
Since all edges are independent, we may write the likelihood over the entire graph by taking a product over all possible edges to give
But the second product is independent of *q* since all possible edges appear and are raised to the same power. Thus, the likelihood of these observations, given any *q*, has the proportionality

So we should choose a reordering *q* to maximize ℒ(*q*). From equation (3.1), we thus have
3.2
Maximizing in equation (3.2) over all reorderings *q* is a hard combinatoric optimization problem. We can make progress through two types of simplification.

First, we assume that the ratio of birth and death transition probabilities has the functional form
3.3
for some constant *θ*. Then taking logarithms in equation (3.2), we have
3.4

The second step is to relax the problem so that *q* is allowed to be a real-valued vector. The right-hand side in equation (3.4) may be written as the quadratic form *q*^{T}Δ_{R}*q*. Here Δ_{R}, the Laplacian matrix associated with *R*, has the form *D*–*R*, where the diagonal matrix *D* contains the row/column sums of *R*. To remove shifting and scaling redundancies we also impose the constraints ‖*q*‖_{2}=1 and .

We now have a tractable optimization problem. For applications where the death rate exceeds the birth rate at long range, so *θ*<1, it is solved by a *Fielder vector*—an eigenvector corresponding to the smallest non-zero eigenvalue of Δ_{R}. A reordering of the vertices can be recovered by sorting the components of *q*; that is, vertex *i* is placed before vertex *j* if *q*_{i}<*q*_{j}. This approach has been found to be effective for static networks (Van Driessche & Roose 1995; George & Pothen 1997; Ding *et al.* 2001; Higham 2003, 2005; Strang 2008; Grindrod *et al.* submitted); here we are showing that in the evolving case it is possible to justify from first principles the idea of reordering on the cumulative (non-binary) matrix *R*.

We point out that it is not necessary to know the actual value of *θ* in equation (3.3). The derivation assumes only that this functional form exists and *θ* is not required by the algorithm. We have found in practice that performance is not sensitive to the precise form of range dependency (Higham 2003), especially in the long range, or large *k*, regime. Furthermore, this approach of spectral reordering based on *R* can be used for any dataset, and the validity of equation (3.3) may then be tested *a posteriori*.

We also note that the reordering approach continues to make sense when *θ*>1 in equation (3.3). This includes the case where long ranges are very unlikely to emerge, but those that do are long lived. In this case, because log(*θ*)>0, the expression in equation (3.4) is maximized by an eigenvector of the Laplacian that corresponds to a dominant eigenvalue.

## 4. Computational results for reordering an evolving graph

To test the reordering approach, we generated some synthetic data from the appropriate underlying model. With *n*=100 vertices, we chose a birth rate *f*_{α}(*k*)=0.1(0.98)^{k2} and death rate *f*_{ω}(*k*)=0.2. After ordering the nodes arbitrarily and letting the network evolve for 200 time steps, we applied the reordering algorithm. In this new ordering, figure 1 shows the binarized sum of adjacency matrices over all 200 time steps; a light dot denotes that an edge was present for at least one time step. Figure 2 shows a typical member of the sequence in this new ordering. We see from figures 1 and 2 that the hidden range dependency has been revealed by the algorithm.

Next, we illustrate the algorithm in an electroencephalography application, using data from Sweeney-Reed & Nasuto (in press). Here, the measurements represent electrical activity produced by the firing of neurons within the brain over a short period of time, reflecting correlated synaptic activity caused by post-synaptic potentials of nearby cortical neurons.

In these experiments, the subject carried out a specific task—tapping in time to music. We use 4 s worth of data sampled at 500 Hz, with 2 s prior to tap and 2 s after. Hence, the finger tap starts around 1000 samples into the data. Measurements were taken at each of 128 electrodes arranged at fixed points on the scalp. Figure 3 shows the electrode locations (see Sweeney-Reed & Nasuto in press for further details).

Results from nearby electrodes may possibly be correlated due to the volume integrative effects of the skin and scalp. On the other hand, transient correlations between channels corresponding to electrodes located some distance apart on the scalp may represent dynamic, synchronous, locking as some otherwise separate tasks become briefly coordinated; or the need for specific tasks that rely on a distributed processing effort. The high-resolution data allow us to examine whether such wide-scale transient correlations occur, as some separate cognitive processes become (briefly) coordinated to manage a combination of processing, sensory and motor responses that may be distributed across the cortex.

We therefore let each electrode represent a vertex within an evolving graph. We subdivided the time-dependent data into windows of 50 consecutive time steps, each lasting 0.08 s. Within each time step we first obtained the *all versus all* channel correlation matrix and defined an edge between vertex *i* and vertex *j* if this correlation exceeded 0.8. This resulted in an evolving sequence of 50 adjacency matrices.

We note that correlation between signals is a far from perfect measure, especially in the search for noisy, transient, synchronous components within time series, but it will serve for our current purpose of illustrating how naturally the concept and methods of evolving (organizational) graphs, introduced in this paper, may represent the coordinated emergent, transient, responses (both local and non-local) of the brain.

Figure 4 shows the evolving sequence of 50 networks. Here time increases along each row in the picture—as shown by the indicative labelling of times 1,2,3,4,5 and 47,47,48,49,50. The effect of the subject’s ‘task’ after 2 s can be seen clearly in time steps 26,27,28. In this figure, vertices are ordered according to the default values provided by the recording equipment, as shown in figure 3, and we note that some vertices with neighbouring indices (i.e. successive row/column indices in the adjacency matrix) correspond to electrodes that are geographically close on the scalp, but in other cases they do not. In particular, there is an artificial periodicity or ‘wrap around effect’ in this default vertex ordering; the earliest vertices are geographically close to the latest. This effect manifests itself through the significant non-zero blocks in the off-diagonal corners of the adjacency matrices.

Figure 5 repeats the information in figure 4, with vertices reordered via the spectral algorithm. In this case, we see that the activity has been arranged into coherent blocks. At the latter end of this new ordering (lower right corner), one set of vertices appears to have a consistently strong set of mutual correlations, whereas at the start (upper left corner) a more transient set of correlations is captured. The apparent periodicity from figure 4 has been removed and there is a clear propensity for ‘short-range’ edges, that is, connections between near neighbours.

In figure 6, we show the binarized cumulative matrix *R* under the new ordering. We also give the new electrode ordering; so electrodes numbered 55, 22 and 17 in the equipment are placed in positions 1, 2 and 3, respectively, and so on. Of particular interest are electrodes 81, 102, 108, 48, 101, 57, 56, 95, 50, 63, 100 and 51, which appear to provide crosstalk between the two otherwise isolated groups. Looking at their physical locations, which are highlighted in figure 3, we see that these regions cover two physically separate areas close behind both ears, and hence it is plausible that these edges correspond to auditory activity.

As an *a posteriori* check on the relevance of the evolving network model, in figure 7 we scatter plot the values
where and are computed from equations (2.2) and (2.3). If assumption (3.3) is valid, then these values provide estimates for . For each predicted range, *k*, the solid line in the figure shows the average of the scattered points, and we see that the results are consistent with the *θ*<1 scenario to which the algorithm applies.

Although it would of course be possible to customize this computational technique to account for the specific nature of the problem, we believe that this experiment on real neuroscience data confirms that the basic modelling and algorithmic approach has potential for understanding, calibrating and summarizing evolving networks. In particular, it allows the transient nature of the interactions to be quantified and compared across different settings.

## 5. Simulating propagation within an evolving graph

In this section, we present two sets of computational results that arise when we put dynamics on the evolving network model. These results are then explained analytically in §6.

### (a) An epidemic model

Suppose we have a binary state variable defined at each vertex and at each time step. To be specific, let this variable take values ‘infected’ and ‘susceptible’. Initially all vertices except one are labelled susceptible.

From one time step to the next, we impose the following simple dynamics, depending on a single parameter *μ*.

A susceptible vertex has no effect on the fate of any other vertex.

Each infected vertex passes on the infection to all of its current immediate neighbours.

Having passed on the infection, each infected vertex becomes susceptible with probability 1−

*μ*, or else remains infected, with probability*μ*.

Intuitively, we foresee two possibilities depending on the infection dynamics and the evolving network dynamics: (i) that the evolving graph is sparse so the infection spreads only for large *μ* (where the infection tends to remain at individual vertices for a long time and these vertices eventually acquire edges), but not for small *μ*, or (ii) the evolving network is dense enough to support spreading for all values of *μ*—even when *μ*=0 the infection lasts for just one time step at each infected vertex. In the special case where *μ*=1 and every *α* is positive, we are in the *successive percolation* or *flooding* (Clementi *et al.*2008, 2009) regime where all vertices are certain to become infected.

The novelty in this area lies in the dynamic coupling between the evolution of the contact network and the time course of the infection, in contrast to most of the existing work in this field, which has been carried out with percolation type models or susceptible, infected, recovered (SIR) dynamics on static graphs.

In figure 8, we depict a threshold case, where the elements of the evolving (range-dependent) network are individually very sparse; the expected degree is 0.79. All elements are highly disconnected. In this example, we use the same evolving network in each case with 100 vertices, for which *f*_{α}(*k*)=0.1(0.9)^{k2} and *f*_{ω}(*k*)=0.5. We seed the infection at a single vertex (vertex 40), with a different value of *μ* for each run. We emphasize that recovery at each infected vertex is decided independently at random in all cases (a Poisson process). We employ values for *μ* between 0.1 (where the infection dies out quickly) and 0.625, where it propagates over the whole evolution. The threshold was observed to be approximately between 0.55 and 0.575. So, we have an example of how occasional long-range transient connections between otherwise short-range, isolated communities may propagate disease providing that the infectious period for individuals (and hence of the isolated communities) is long term enough.

### (b) A message passing viewpoint

Suppose that *μ*=1, and that we are transmitting a message rather than an infection. Starting out from a specific vertex, called the sender, at each time step the message is relayed to any vertices directly connected to previous message holders. Hence, the message spreads according to the successive edges added at each time step to the existing message holders. Eventually, all vertices will receive the message (since *μ*=1 so a vertex never forgets the message). Suppose also that a vertex is designated as the desired *receiver*. Then, in the range-dependent graph setting, there is a natural range between sender and receiver, which we will denote by *k*_{mess}. We will say the minimum number of time steps needed for the first arrival of the message at the receiver is the *shortest transmission time* (STT). It is clearly of interest to study the distribution of the STT and related quantities, such as the natural ‘path length’ measured by taking whichever message-holding vertex is currently closest in range to the receiver. This path length is important if, for example, the message propagation from one vertex to another along any single edge during any single time step is noisy, so that a ‘Chinese Whisper’ effect accentuates such noise.

For *f*_{α}(*k*)=0.1(0.9)^{k2} and *f*_{ω}(*k*)=0.5, we have an average degree at any time step given by *z*=0.792. For *n*=200 we selected both a sender vertex and a receiver vertex randomly and carried out an experiment with a fresh evolving graph in each case. Then over 17 300 such experiments, we obtain the results shown in figure 9 for the mean, mode and distribution of the STT.

## 6. Analysing propagation within an evolving graph

We now show that the evolving random network model used in figures 8 and 9 is sufficiently compact to allow for some theoretical analysis.

### (a) Analysis of the message passage simulations

In order to understand the behaviour seen in figure 9, we will make two simplifying assumptions.

At each time step, the nearest vertex to the receiver, say

*A*, either remains the nearest vertex, or the new nearest vertex arises from an edge that has appeared from vertex*A*. We may then assume that each new edge utilized is independent of the previous ones (since we make no assumption about any edge leaving*A*, prior to the arrival of the message at*A*).The ‘edge effects’ caused by the requirement that the message exactly reaches the receiver can be ignored, so that we can focus on the general phase where the message progresses as quickly as possible away from the sender.

In any element of the evolving graph, we have the equilibrium probability that any edge of length *k* is present at a particular time step given by equation (3.1).

Consider any particular vertex that the message has reached, and let *π*(*k*) denote the probability that the longest range of any edge connected to it is exactly equal to *k*. Then, we have directly

The sequence {*π*(*k*)} sums to one, and we can use these forms to calculate the expected range of the longest range edge connected to any vertex .

Similarly, let *π*^{★}(*k*) denote the probability that the longest range of any edge connected to the particular vertex *in the direction of the receiver* is exactly equal to *k*. Then, we have

The expected range of the longest range edge (LRE) in the direction of the receiver is thus .

For the choices of *f*_{α}(*k*) and *f*_{ω}(*k*) in figure 9, the expected range of the longest range edge is given by 1.3598 in any direction. Similarly, the expected range of the longest range edge in the direction of the receiver is given by LRE=0.773.

Under our simplifying assumption that edges transmit the message in the direction of the receiver, we have

This explains the slope of the mean and mode curves in figure 9. Over long message ranges we may appeal to the central limit theorem: the message propagates according to a process that could be subdivided into independent increments, for much shorter message ranges, successively drawn from a distribution of paths for the shorter ranges. Thus, the message arrival behaviour and time resembles that for a diffusion–advection process.

Of those successive time steps, we expect that a fraction *π*^{★}(0)=0.657 made no progress towards the receiver—the range of the longest edge in the receiver’s direction being zero—for that time step at the closest vertex.

Hence, the message travels a total message range of *k*_{mess}, over 1.293 *k*_{mess} successive time steps, along approximately 1.293(1−0.657)*k*_{mess}=0.442*k*_{mess} edges. Notice that the expected range of the longest range edge in the direction of the receiver, given that it is not zero, is simply , which is approximately 2.25.

### (b) Analysis of the epidemic simulations

For the computations shown in figure 8, we observed a threshold value for *μ* governing overall growth or decay of the disease. To understand this behaviour and estimate the critical *μ*, we consider, at equilibrium, a single infected vertex, *A*. Our aim is to estimate the number of new vertices that *A* infects before its eventual recovery. The expected degree of node *A* on the next step is the sum over the probability of all possible edges,
6.1
where *p*_{0}(*k*) is defined in equation (3.1). Some of the vertices that *A* points to may already be infected, of course. In particular, the most likely already-infected neighbour is the vertex, say *B*, that infected *A*. Because the edge death rate, *f*_{ω}, is constant in this experiment, this value gives the probability that the edge connecting *B* to *A* remains over one step. So the probability that the edge *A*–*B* remains *and* vertex *B* has not recovered is *f*_{ω}*μ*. Hence, correcting in equation (6.1) to discount ‘infection’ of an already-infected node *B*, we arrive at
as our approximation to the expected number of new infections caused by *A* after one step.

Now let us consider new edges appearing at subsequent steps. Suppose we consider any such edge of length *k*. This edge is not present on the first time-step post-infection with probability 1−*p*_{0}(*k*). The probability that it *first* appears as a new edge on exactly the *j*th time-step post-infection (*j*=2,3,…) is thus
(This expression arises as the product of the probabilities that the edge (i) is missing on the first step, (ii) does not appear in steps 2 through to *j*−1, and (iii) appears at step *j*.) The expected number of edges connecting to vertex *A* for the first time at the *j*th time-step post-infection is thus

The probability that vertex *A* is still infectious at the *j*th step following infection is *μ*^{j−1}. So, summing over all times, our overall approximation of the expected number of new vertices that *A* infects before recovery is
6.2
This quantity is analogous to the basic *reproduction rate*, usually called *R*_{0}, that arises in classical epidemiology (Diekmann *et al.* 1990), and we are therefore interested in the critical value *R*(*μ*)=1.

For the case *f*_{ω}=0.5, *f*_{α}(*k*)=0.1(0.9)^{k2} used in figure 8, we show in figure 10 the behaviour of *R*(*μ*) in equation (6.2). Solving numerically, we find that *R*(*μ*)=1 at *μ*=0.575, which matches well with the experimentally observed threshold.

## 7. Discussion

Many new challenges arise when we move from fixed networks to the more general case where connections are time dependent. In this work, we highlighted some fundamental tasks in the modelling and calibration of evolving networks and proposed a novel range-dependent birth and death mechanism. This framework permits a variety of evolutionary behaviours but also imposes enough structure to make it possible for some analysis and algorithm development. We showed that the evolving network model adds value to some real brain activity data, but we believe that, due to the generic nature of the modelling assumptions, the same approach will be applicable to many other scenarios where sparse network ‘snapshots’ are observed, including applications in telecommunications, sociology and business.

Although our evolving network model was set up as a discrete-time Markov chain, we remark that an analogous continuous time framework could be developed, which might be more suited to rigorous probabilistic analysis.

The evolving network model can also be combined with a time-dependent process that acts *over* the network. We showed via simulation and analysis that the overall behaviour of the resulting stochastic process depends strongly on the interaction between the two sets of dynamics. There is much scope here for adding realistic topological dynamics to the traditional static network models of disease propagation and message passing.

## Acknowledgements

We thank Slawomir Nasuto for making available some data from Sweeney-Reed & Nasuto (in press) and a template for figure 3, and for helpful discussions about the computational results in §4. The authors acknowledge support from EPSRC through project grant GR/S62383/01 and Bridging the Gaps ‘Cognitive System Science’ grant EP/F033036/1, and from MRC through project grant G0601353.

## Footnotes

- Received August 27, 2009.
- Accepted October 13, 2009.

- © 2009 The Royal Society