## Abstract

In this paper a random link activation–deletion (RLAD) model is proposed that gives rise to a stochastically evolving network. This dynamic network is then coupled to a simple susceptible-infectious-suceptible (*SIS*) dynamics on the network, and the resulting spectrum of model behaviour is explored via simulation and a novel pairwise model for dynamic networks. First, the dynamic network model is systematically analysed by considering link-type independent and dependent network dynamics coupled with globally constrained link creation. This is done rigorously with some analytical results and we highlight where such analysis can be performed and how these simpler models provide a benchmark to test and validate full simulations. The pairwise model is used to study the interplay between *SIS*-type dynamics on the network and link-type-dependent activation–deletion. Assumptions of the pairwise model are identified and their implications interpreted in a way that complements our current understanding. Furthermore, we also discuss how the strong assumptions of the closure relations can lead to disagreement between the simulation and pairwise model. Unlike on a static network, the resulting spectrum of behaviour is more complex with the prevalence of infections exhibiting not only a single steady state, but also bistability and oscillations.

## 1. Introduction

Many real-world systems ranging from neuroscience and epidemiology to computer sciences and socioeconomics can be represented as well-defined units interacting via a static or dynamic set of links or connections (Broder et al. 2000; Kossinets & Watts 2006; Vernon & Keeling 2009; Sporns 2011). The wide applicability of networks as a modelling tool has captured the attention of many different research communities and has led to the development of a large body of research at the interface of network/graph theory, stochastic processes, probability theory, discrete mathematics and computer sciences (Albert & Barabási 2002; Newman 2003). Initially, most of the research concentrated on the structure and properties of real-world networks and aimed to understand and uncover the laws that gave rise to the observed or empirical networks. In parallel, a different research direction emerged, namely the study of how the properties and structure of the network impact on the dynamical processes taking place on it (e.g. flow of information on the Web, disease transmission on social networks and self organization of neurons). While earlier research focused on the dynamics or evolution of networks (Dorogstev & Mendes 2002; Chan et al. 2003) without considering the dynamical processes they support, latter research considered fixed and static networks and mainly focused on the dynamics on networks. However, in many cases considering both the dynamics of the network and on the network is essential to understand the problem that is being modelled, but doing so raises several challenges. Firstly, capturing and modelling the interaction between the two dynamics are non-trivial with little empirical evidence and second, the increased complexity reduces analytical tractability with results mainly relying on simulation.

Over the last few years, various simple classes of models have been proposed where both the dynamics of the network and on the network are considered (Gross & Blasius 2008). For example, Saramäki & Kaski (2005) proposed a model where nodes are distributed on a ring and links are divided into short-range (SR) and long-range (LR) links. SR links are considered to be fixed and connect nodes to their nearest neighbours. LR links vary randomly, meaning that an infected node tries to infect with some probability, a node chosen at random and succeeds to do so if the chosen individual is susceptible. The authors formulated a simple ODE model that is similar to pairwise models and used this to derive analytical results for disease transmission threshold and to validate simulation results.

Gross et al. (2006) proposed a network-based model where (*SI*) links/pairs (i.e. a link between a susceptible (S) and infectious (I) individual) are broken at a certain rate with susceptible nodes immediately re-wiring to other susceptible nodes chosen at random from the entire population. Again, a simple pairwise type model was used to derive a low-dimensional dynamical system that describes the interaction of network and disease dynamics. This model operates on the strong assumption that the status of each node is globally available. Risau-Gusman & Zanette (2009) relaxed the assumptions of the model proposed by Thilo *et al.* and considered the case when the susceptible node from an (*SI*) pair re-wires to a node chosen at random from the entire network regardless of its state. They generalized the model further and considered the case where not only susceptible nodes are allowed to change their contacts, but also the infectious nodes can do so. In this model, when an *SI* link is cut, the *I* node will reconnect with probability *q* to a node chosen at random and independently of its state, and with probability 1−*q* the susceptible node will keep the contact and look to reconnect at random. Both papers use pair-approximation models to validate simulation results.

Grindrod & Higham (2010) proposed two simple models of undirected evolving graphs by starting from the complete state space with possible elements or different graphs over *N* nodes. With this in place, the evolution of a graph can be represented as a path in the state space governed by some stochastic process *P*(*G*_{i+1}|*G*_{i},*G*_{i−1},…), where *G*_{i} is the state of the graph, represented as a symmetric adjacency matrix, at discrete time step *i*. While this formulation is exact and can be extended to continuous time, the drawback of this approach comes from the large number of transition probabilities for time independent processes (i.e. 2^{N(N−1)}) or the Kolmogorov transition equations (i.e. 2^{N(N−1)/2}). The authors proposed an edge birth-and-death model where each edge can become activated or deleted at a fixed probability that is independent of all other edges or the current state of the graph. They extended this to a more sophisticated model where link activation and deletion depend on some form of proximity between edges as given by the initial set-up of the nodes (e.g. fixed location on a line with nodes connected according to some connectivity kernel). The aim of their study was twofold: first to use such models to simulate a series of dynamic or evolving networks and second to fit evolving network models to data and use likelihood-based estimation of model parameters.

Miller et al. (2011) and Volz & Meyers (2007), using the probability generating function (PGF) formalism or edge-based compartmental modelling, have recently also studied a variety of dynamic network models (e.g. mean-field social heterogeneity (MFSH) and dynamic fixed degree with or without dormant contacts) coupled with susceptible-infectious-recovered (*SIR*) dynamics. However, these network dynamics are not coupled with the epidemic and the authors propose low-dimensional ODEs that capture the time evolution of the epidemic well. Calculations of the basic reproductive number, *R*_{0}, are also presented with focus on the impact of partnership duration and social heterogeneity on epidemics.

In this paper, we consider a set of simple dynamic networks based on the link birth-and-death or activation–deletion model, and also a variant of this, where link birth/activation is globally constrained. In order to provide a succinct overview of different analytical methodologies and their limits when applied to such models, the paper considers four scenarios of increasing complexity: the impact of simple network dynamics on the structure of the network when node dynamic is absent (§2, unconstrained and globally constrained link creation), when the nodes are static but labelled (§3) and when the dynamic network is coupled with susceptible-infectious-suceptible (*SIS*) node dynamics (§4). To conclude, we study in depth the impact of network dynamics and its coupling with the disease transmission on the outcome of epidemics (§5). For the first three models, results are obtained by using the exact formulation in terms of Markov Chains and Kolmogorov equations and a mean-field type approach with further results from simulation. These methodologies cannot deal with the final scenario, however, the output of the analytical models can be used as a benchmark to test the validity of the full simulation model. The full simulation model is used in conjunction with a novel pairwise approximation to investigate and characterize the full spectrum of behaviour. The agreement between simulation and pairwise model is discussed in detail. This is done by providing some new insight into how the assumptions of the pairwise model can be interpreted and how the strong assumptions of the closure relation can lead to disagreement between the two models. Parameter regimes where agreement is excellent are identified, however, many open questions remain regarding the exact relation between simulation results and output from the pairwise model.

## 2. Simple models of stochastically evolving networks

In this section, the impact of simple network dynamics on the structure and properties of the network is explored. We present different modelling techniques to derive exact or approximation models and compare these to simulation results. Two simple network dynamics are considered: (i) random link activation–deletion (RLAD), whereby non-active or non-existing links are activated with a given rate (e.g. *α*), while existing ones are deleted with some other rate (e.g. *ω*) and (ii) globally constrained RLAD, whereby the above link creation process is constrained globally, i.e. the higher the number of active links the lower the rate at which new links are activated. Initially, the simple *dynamics of the network* are studied without considering any *dynamics on the network* or any link activation or deletion that may depend on node labels.

### (a) Random link activation–deletion: probabilistic approach

Let us consider an undirected and unweighted network with *N* nodes where the maximum number of edges is *M*=*N*(*N*−1)/2. The dynamics of edges evolve according to the following two simple rules. Non-active or absent links are activated independently at random at rate *α* while existing links are broken independently at random at rate *ω*. Let (*X*(*t*))_{t≥0} be an integer-valued random variable that represents the number of edges/links in the network at time *t*. If *P*(*X*(*t*)=*k*)=*p*_{k}(*t*) then the Kolmogorov equations for *p*_{k} are given by
2.1where *k*=0,1,2,…,*M* with *p*_{M+1}(*t*)=*p*_{−1}(*t*)=0 (no link activation when *k*=*M* and no link deletion when *k*=0). The average number of edges in the graph at time *t* can be defined as . Upon using equation (2.1) and after some simple algebra, the equation for *K*_{1}(*t*) follows easily and is given by
2.2Hence, at equilibrium
The equilibrium value can also be found heuristically by determining the number of edges at which the total rate of link activation () is balanced by the total rate of link deletion (). For this set-up, the average number of links per node or average node degree is given by

A traditional and commonly used method to determine the probability distribution as given by the forward Kolmogorov equations above (equation (2.1)) is via the PGF technique. This is defined as
2.3where *s* is a placeholder variable that allows to concentrate all the information about the probability distribution into one single function. Multiplying equation (2.1) by *s*^{k}, followed by summation for *k*=0,1,…,*M* and some simple calculations gives rise to the following partial differential equation
2.4and
2.5where the initial condition corresponds to starting with *m*_{0} edges (i.e. *p*_{m0}(0)=1 and *p*_{k}(0)=0 for ∀*k*∈{0,1,2,…,*M*}\{*m*_{0}}). This is a first-order homogeneous partial differential equation that can be solved using the method of characteristics and its solution is given by
2.6From the equation above, *K*_{1}(*t*), the expected number of edges in the network is given by , and in the limit of , *K*_{1}(*t*) tends to its equilibrium value defined as
2.7Moreover, for
2.8which is the generating function of the binomial distribution with *M* trials and a per-trial probability of success *p*=*α*/(*α*+*ω*). The excellent agreement between the exact model and simulation is shown in figure 1*a*, and confirms that the formulation of the two models is correct and consistent. The number of equations in the exact model scale as *O*(*N*^{2}) which restricts its applicability to small networks.

### (b) Random link activation–deletion: mean-field approach

To overcome the limitations of the exact model and to gain more insight into the properties or structure of the network, a different modelling approach is presented below. Here, the network is considered as a population of nodes where nodes can be classified according to their number of links/contacts and the rates of moving to either more highly or less well connected classes. Hence, the modelling relies on deriving evolution equations for the number of nodes with degree 0,1,2,…,*N*−1. Let *n*_{k} denote the number of nodes with *k* contacts where 0,1,2,…,*N*−1. Based on simple heuristic reasoning, the evolution equations for *n*_{k}s are given by
2.9This approach is similar to compartmental models where, in this case, transitions between compartments represent link gain or link loss. Solving the equations above provides information about the number of nodes with different connectivity but does not test whether such a network is realizable. This bears strong similarities to some modelling approaches where networks are simply considered in terms of nodes and their stubs but without being connected up in a coherent network (Ball & Neal 2008; Lindquist et al. 2011). The linear system of ordinary differential equations lends itself to some simple analysis. The rates matrix is of particular interest as this determines the eigenvalues (*λ*_{j}) and the corresponding eigenvectors (*s*_{j}), where *j*=0,1,…*N*−1, that are used to construct the solution of the system. The solution of the system in terms of these is given by
where *c*_{j}s are some arbitrary constants. Owing to the special structure of the rates matrix, all eigenvalues but one have negative real parts. The single eigenvalue with a *non-negative real part* is in fact a zero eigenvalue. Hence, when , the solution of equation (2.9), *n*(*t*)=(*n*_{0}(*t*),*n*_{1}(*t*),…,*n*_{N−1}(*t*)) tends to the eigenvector associated with the 0 eigenvalue times some constant, when . It is easy to show analytically that
This is equivalent to the binomial distribution with *N*−1 trials and a per trial probability of success given by *α*/(*α*+*ω*). This is further illustrated in figure 1*b* where the degree distribution, for large time, is plotted for different 〈*k*〉 values. Again, the agreement of the mean-field model with simulation is excellent, and confirms that, as expected, the simple RLAD model leads to simple Erdos–Rényi type networks.

### (c) Globally constrained random link activation–deletion: probabilistic approach

Both real and theoretical networks are usually sparsely connected with 〈*k*〉≪*N* and thus it is natural to impose a limit on the total number of edges in the graph. To do this, the simple RLAD model is extended to include a carrying capacity () that can limit link activation. This carrying capacity can be viewed as the maximum number of links allowable in the network. Of course this is not the only approach to limiting the number of edges in the network. It could be argued that a more realistic approach would be to introduce a local neighbourhood capacity to nodes in the network (see Taylor et al. 2012). Here, we continue with the global constraint and following the same notation as above, the globally constrained RLAD (GC-RLAD) model leads to the following equations
2.10where *k*=0,1,2,…,*M* with obvious modifications for *k*=0 and *k*=*M*. Using a similar approach as above, the equation for *K*_{1}(*t*), the average number of edges in the graph at time *t*, is given by
2.11where is the second moment of the edge number distribution at time *t*. This equation cannot be solved directly since it involves the second moment of the edge distribution and an equation for this is needed. As an alternative to this exact approach, we consider the most obvious and simplest approximation whereby . The following theorem shows that this approximation becomes exact in the limit of .

### Theorem 1

*For the solution* *of equation (2.11) upon substituting K*_{2}*(t) by* *the following statement holds. Provided that* *remains constant as* *then for any T>0 there exists a constant C>0 such that
**where* *.*

The proof of the theorem is rather technical and it will be presented in a more wider context as part of an additional paper. However, figure 2*a* confirms our statement and shows the excellent agreement between the exact and the approximation model even for small networks, where *M* is also small. Using this result, a quadratic equation for the edge value at equilibrium () follows easily,
2.12While both solutions are positive, the smaller of the two can be confirmed as the correct equilibrium value given that the number of edges in the network cannot be larger than the carrying capacity. In the case of , our argument follows immediately.

### (d) Globally-constrained RLAD: mean-field approach

Following the same considerations as for the RLAD model and using the same notation, the globally constrained RLAD at the node level is characterized by the following equations
2.13where *e*(*t*) is given by
In figure 2*b*, for a given set of parameters, we show that the mean-field approach can again be used to predict the degree distribution of the network, and as before, the numerical results suggest that networks based on the GC-RLAD dynamics are also randomly connected with binomial distribution. However, in this case, the good agreement with the binomial distribution could not be demonstrated analytically directly from equation (2.13).

## 3. Simple stochastically evolving networks in the presence of node labelling

The network dynamics so far has only been considered without the dynamics on the network and independently of node and/or link type. To explore the full impact of the dynamics on the network, it is advisable that a step-by-step approach is taken as proposed by Gross et al. (2006). First, the impact of network dynamics can be investigated independently of node and/or link characteristics which could either be imposed externally and be static or be dynamic as a result of a separate dynamics running on the network. Second, the network dynamics based on static node/edge labelling can be investigated, and finally, the impact of coupling the dynamics of the network with the dynamics on the network needs to be considered. In the latter and most interesting case, the state of the nodes will have an impact on the network dynamics and vice versa (e.g. breaking susceptible-infected links impacts on the networks with the converse also being true; neuronal activation impacts on networks through experience-dependent plasticity and vice versa).

### (a) SI labelling

Be it disease or information transmission, the manipulation of networks via preferential node/link addition and/or deletion provides a powerful mechanism to influence and optimize processes unfolding on the network. For example, disease transmission can be slowed or halted if links between susceptible and infected individuals are cut fast enough (Gross et al. 2006). Using the analogy of simple epidemic models, such as the *SIS* model, nodes are labelled at random as *S* or *I*. In the spirit of the network dynamics so far, all pair types can be activated or deleted at random with pair-type-dependent rates. The rates at which (*SS*), (*SI*) and (*II*) links are activated are denoted by *α*_{SS}, *α*_{SI} and *α*_{II}, respectively. Similarly, *ω*_{SS}, *ω*_{SI} and *ω*_{II} represent the rates at which (*SS*), (*SI*) and (*II*) links are deleted. These labels are permanent and do not change during the evolution of the network. In this case, the Kolmogorov equations can be written down but are far more complicated compared with the previous cases. The state space is now given in terms of three variables, namely the counts of the various link types (*SS*), (*SI*) and (*II*). These take values from 0 to {*S*}({*S*}−1), {*S*}{*I*} and {*I*}({*I*}−1), respectively, where {*S*} and {*I*}=*N*−{*S*} are the initial number of nodes labelled *S* and *I*. It is more practical to use the mean-field approach and write down ODEs for the number of different pair types. For brevity, only the equations for the globally constrained case are given,
3.1
3.2
and
3.3where now *K* is the carrying capacity, and *e*(*t*) is given by
and represents the number of doubly counted edges. It is worth noting that here *n*_{AB} (where *A*,*B*∈{*S*,*I*}) stands for doubly counted edges (i.e. a single *SS* pair counts for two *SS* edges in *n*_{SS}). Similarly, the carrying capacity *K* needs to be understood as such. The expected number of pairs at equilibrium (*SS*_{eq},*SI*_{eq},*II*_{eq}) are given as the solutions of the following simple equations,
3.4
3.5
and
3.6where
and {*S*}({*S*}−1), {*S*}{*I*} and {*I*}({*I*}−1) are constants determined by the initial number of *S* and *I*s. These equations allow us to find the equilibrium values at which the total rate of link activation equals the total rate of link deletion. The following proposition shows that the solution of equations (3.4)–(3.6) is unique.

### Proposition 1

*The system given by equations* (*3.4*)–(*3.6*) *has a unique solution given as*
*where*
*and*
*and y is the unique root of the function*
*in the interval* **[0, K]**.

*Finding the solution of the equation h*(

*y*)=0

*reduces easily to find the root of a fourth degree polynomial*.

### Proof.

Multiplying equations (3.4)–(3.6) by *K*/*α*_{SS}, *K*/*α*_{SI} and *K*/*α*_{II}, respectively, and introducing the following new notation
allows us to recast equations (3.4)–(3.6) as follows,
Substituting these expressions into the definition of *y* and using that
the following equation for *y* is obtained
This shows that *y* is a root of function *h*. It is easy to see that *h* is a decreasing function, hence it has at most one root. Moreover, *h*(0)=*K*>0 and *h*(*K*)<0, since *K*_{i}*A*_{i}/(*A*_{i}+*K*)<*K*_{i}. Therefore *h* has a unique root in the interval [0,*K*].

It is important to note that if the network dynamics is not globally constrained, equations (3.4)–(3.6) decouple and each equilibrium value can be found based on results from the non-constrained random link creation–deletion case. The excellent agreement between the theoretical equilibrium, as defined above, and simulation is illustrated in figure 3.

## 4. Interaction of network and disease dynamics: simulation and pairwise model

In this section, the dynamics of the network and the simple *SIS* model are studied concurrently using individual-based stochastic simulations and a low-dimensional, deterministic pairwise model that acts as an approximate model in the sense that, under appropriate conditions, output from the stochastic/exact model should tend to results based on the deterministic limit.

### (a) Simulation model and validity checks

The simulation is implemented on a network of fixed size with *N* nodes. Disease transmission (*τ*-transmission rate across an (*SI*) link), recovery (at rate *γ*), link creation (at rate *α*_{ab} with *a*,*b*∈{*S*,*I*}) and deletion (at rate *ω*_{cd} with *c*,*d*∈{*S*,*I*}) all act simultaneously and give rise to the coupled disease and network dynamics. The simulation is based on synchronous updating, where, to be in line with the Markovian framework, the simulation time step (Δ*t*) is kept small enough to try and ensure that not more than one event/update happens per single iteration. In a small time interval Δ*t*, a susceptible node with *k* infectious neighbours becomes infected with probability , while an infectious node recovers with probability . The activation and cutting of links can be easily performed by considering all pairs of nodes, independently of whether they are connected. If connected, an edge of type (*ab*) is cut with probability , while a non-existing edge between nodes of type *c* and *d* is activated with probability . For globally constrained link creation, the unconstrained rate (*α*_{cd}) is simply multiplied by (1−*f*), where *f* is the fraction of existing over allowable links given by the carrying capacity. All the above processes are modelled as simple Poisson processes that happen independently of each other.

The resulting simulation has an added level of complexity compared with simple disease transmission models, and it is therefore important that the simulation model be validated. In this case, the checking and benchmarking of the simulation are possible owing to the theoretical work developed in §§2 and 3. For example, to test the correctness of the implementation of the link-activation deletion mechanism, the following checks have been performed and completed successfully: (a) comparison of simulation output, with no epidemic and link-type independent activation–deletion, to analytical results from §2, (b) comparison of simulation output, with no epidemics but with labelled nodes, to analytical results from §3. In both cases, the basis for comparison was the average number of links, and in the case of (b), the average number of different link types. While in this paper, the globally constrained RLAD was not coupled with epidemics, it was implemented and tested against analytical results from §2.

### (b) Pairwise model formulation for a dynamic network

In line with the standard pairwise models (Keeling 1999; Simon et al. 2011) and the notations established above, we can now heuristically write down equations for the rate of change of individuals and pairs for the dynamic network as follows,
4.1
4.2
4.3
and
4.4The system of equations is complemented by theoretical initial conditions or taken as expected counts from the simulation model. Triples are closed according to the simple closure given by
where *n* is the average number of links per node. This closure implicitly assumes that the precise type of nodes (e.g. *S* or *I*) around a node in state *B* are independent. This type of closure is widely used for homogeneous random graphs (Keeling 1999; House & Keeling 2010) with no clustering and also for unclustered random graphs with close to Poisson degree distribution (Taylor et al. 2011). The equations above account for link activation–deletion and, without any further constraints, have been made dependent on the type of link being activated/deleted. For example, terms such as *α*_{SI}((*N*−[*I*])[*I*]−[*SI*]) and *ω*_{SI}[*SI*] account for the activation and deletion of (*SI*) links, where (*N*−[*I*])[*I*]−[*SI*] denotes the number of potential (*SI*)-type links that are not yet connected.

### (c) Comparison of pairwise model to simulation: a general consideration

In formulating such a low-dimensional model, the immediate key question is whether such a simple system can approximate, at least qualitatively, results obtained from the simulation model. The answer to this question depends strongly on the network structure, the type of dynamics and how the proposed closure performs when these factors combine. An important observation that is generally valid for all pairwise models with the known closures is that terms such as −*τ*[*SI*] are exact. Thus, the network structure and the formation of correlations as measured by , where *A*,*B*∈{*S*,*I*}, should be conserved. However, the evolution equation for [*SI*] relies on the exact expectation of triples which is approximated by the closure and thus is not exact. We note that values of the correlation measure close to one indicate that the true expected number of pairs is close to the value obtained in the case where nodes are labelled at random as *S* and *I* rather than as a result of the dynamics unfolding on the network. The strong assumption of independence in the closure leads to a closed ODE system that dissipates the true correlations to some degree but, not completely. This argument follows simply from considering a simpler closure at the pair level, i.e. [*SI*] can be approximated by [*SI*]=*n*[*S*][*I*]/*N*. This obviously removes all correlations and gives rise to a system that is equivalent to the mean-field *SIS* model. Using the closure at the level of triples, rather than pairs, will conserve some of the correlations at pair level, however, a deviation from the simulation model is unavoidable. To shed some further light on the impact of closure, it is important to make a distinction between whether nodes are labelled at random as *S* and *I* or whether (*SS*), (*SI*) and (*II*) links are placed at random to form a labelled graph (labelling is to be understood as carried out at random and fixed in time as opposed to being the result of an epidemic unfolding on the graph). Nodes labelled at random will lead to all pair correlation being close to one, , where *A*,*B*∈{*S*,*I*}, in which case, the triples can be counted in terms of singles as follows:
This in turn is equivalent to the closure in terms of pairs. The main difference between the two closures is that the first, more realistic one, allows for pair correlations that are different from one, and thus not random, while keeping the distribution of links random. Thus, the first closure can accommodate more realistic scenarios where non-random pair correlation and random link distribution can coexist. For the agreement between simulation and ODE (equations (4.1)–(4.4)) to hold, we can give the following sufficient conditions:

— , where

*A*,*B*∈{*S*,*I*},— [

*ABC*]=((*n*−1)[*AB*][*BC*])/*n*[*B*] holds and the pair correlation can differ from one, , where*A*,*B*,*C*∈{*S*,*I*},

where the first implies the second but not vice versa.

The conditions above are rather strong and are easily violated even for simple static networks where the ODE model with an even more sophisticated closure (House & Keeling 2010) will fail to match the outcome of simulation as shown by Taylor et al. (2011). However, these conditions are not necessary to get good agreement. Even when the second condition is not fulfilled (i.e. [*ABC*]≠((*n*−1)[*AB*][*BC*])/*n*[*B*]), it is still possible to get reasonably good agreement over time in the expected number of infecteds as illustrated in figure 4.

Our discussion of the results and agreement between the two models will revolve around the arguments presented above and the concept of *preferential link or edge dynamics*. These will be used to underpin our explanation of results from model comparisons. Preferential link dynamics is a direct result of the model ingredients, which allows us to tune link activation and deletion such that certain type of links can be over- or under-represented compared with a random link distribution scenario. Moreover, this means that nodes of certain type preferentially connect to nodes of similar or different nature leading to non-trivial correlation structures on the network. This alone can explain why the ODE with the triple approximation cannot capture such type of correlation structures on the network. In extreme cases, this can be easily illustrated by regimes where the graph breaks down into sub-graphs that are dominated by either *S* or *I* nodes and where these sub-graphs can be close to fully connected or with no edges. For example, if out of all link activation–deletion rates only *ω*_{SS} and *α*_{II} are non-zero, it is straightforward to obtain a regime where the network over time is such that *S* nodes are isolated and lose all their links, while the *I* nodes tightly cluster into an almost completely connected graph. In cases, such as these, the ODE will fail to capture the true dynamics owing to the isolation of susceptibles.

While correlations will not necessarily invalidate the agreement between the two models we identify them as an important factor in determining whether the two models will agree. Other important factors can contribute to whether good agreement is observed. For example, the parameter values could be such that an initially unclustered network may become clustered with higher link density and in this case, the assumptions of the simple closure will break down. The effect of clustering, however, can be counteracted by the network becoming even more densely connected, at which point, the mean-field limit can be approached.

## 5. Model behaviour: impact of network dynamics on epidemic

In this section, the focus is on identifying possible system behaviours and understanding the coupled impact of disease and network dynamics on epidemics. The analysis starts with the consideration of the closed pairwise system which is four dimensional and lends itself to standard bifurcation analysis to determine all possible steady states and their stability as a function of the model parameters. Using the parametric representation method (PRM; Simon et al. 1999), a more rigorous approach used to investigate global bifurcations, four different regimes are identified (figure 5). The analysis shows that the interaction of network and disease dynamics leads to a spectrum of behaviours ranging from a single stable disease-free steady state to stable oscillations. More precisely, there are three types of bifurcations. First, a transcritical bifurcation where the disease-free steady state loses stability with the disease becoming established and thus giving rise to a stable endemic equilibrium. Second, a saddle–node bifurcation which gives rise to the coexistence of two stable equilibria (one being disease-free and the other endemic) with an unstable equilibrium; and finally, a Hopf bifurcation, where the stable endemic equilibrium looses its stability and gives rise to a stable limit cycle. From a disease control viewpoint, of main interest is the region above the transcritical and saddle–node bifurcation curves. In this region, only the disease-free equilibrium is stable indicating that the deletion or breaking of (*SI*) edges curtails the spread of the epidemic and leads to a desirable disease-free steady state. As expected, the regions where different model behaviours are observed can vary and depend on model parameters. For example, in figure 5, the region of bistability is small compared with the Hopf island that covers a considerable part of the parameter space. This case was chosen to illustrate a range of potential outcomes and the rich behaviour of a relatively simple model that can lead to outcomes other than a stable disease-free or stable endemic state, as in the case of the simple *SIS* model on a static network.

Owing to the dynamic nature of the network, the average degree of the nodes (i.e. *n* as given in the triple closure) is a variable itself and changes with time. Hence, the analysis above performed for a fixed *n* serves only as an indicator of possible system behaviours but can give good results if *n* is a slow variable where, for example, the network dynamics is much slower compared with the epidemic or if *n* does not vary considerably. We also note that *n* only enters via the (*n*−1)/*n* term which for realistic networks that are well connected is close to one. For dynamic networks with type-dependent rates of link activation and deletion, the average degree of a susceptible node can be different to that of infective nodes. It is therefore more appropriate to use the following triple closure
where *n*_{B} is the average degree of a node of type *B*. Replacing *n* by *n*_{S}(*t*)=([*SS*]+[*SI*])(*t*)/*S* in the ODE prevents the straightforward derivation of analytical results but allows for fast numerical exploration of the parameter space to identify different possible model outcomes. If in good agreement with simulation, the ODE provides a simple and easy way to explore the full spectrum of behaviour for a large number of parameters, a problem that is extremely difficult to tackle via simulation alone.

The interplay between the two dynamics merits further scrutiny from the viewpoint of the relative time scale of the two processes. In particular, three different regimes can be identified: (i) fast network dynamics compared with the spread of the epidemic, (ii) comparable time scales, and (iii) slow network dynamics relative to disease transmission. As expected, in the limit of networks that mix fast (i.e. regime (i)), the simulation approaches the standard mean-field model, figure 6, given by
where 〈*k*〉 represents the average degree, once the network has reached its equilibrium. For quick network dynamics, 〈*k*〉 is approached much faster compared with the time scale of the epidemic, and this value can be used as an input in the standard mean-field model above. This approach is in fact equivalent to decoupling the system of two ODEs composed of the simple mean-field model above and the evolution equation for *K*_{1}(*t*) given by equation (2.2). For fast network dynamics, 〈*k*〉 quickly approaches 2*K*_{1}/*N*, where with creation and deletion rates chosen to obtain a plausible, positive value. Thus, the average degree of the network at the steady state serves as input to the mean-field epidemic model and the system decouples. The setup in the third regime leads to some interesting behaviour of the coupled system. A slow network dynamics, relative to the spread of the epidemic, may make the coupled process look like disease transmission on a static network. However, this will only hold initially as shown in figure 6, where a well-established epidemic can be brought to a halt by a slow but potent network dynamics that thins the networks sufficiently. The same figure also illustrates that when the thinned network can still support an epidemic, the slow network dynamics applies a delayed but severe correction to the steady state, i.e. reduced infection prevalence level. When the timescales are comparable, the processes evolve hand-in-hand and this regime is explored further both from an application and model agreement viewpoint.

Good agreement between pairwise model and simulation would demonstrate that the ODE model can be used as an effective tool to investigate possible model outcomes. In what follows, we identify parameter sets where such agreement for the coupled dynamics case exists. The results are based on extensive individual-based stochastic network simulations for a large selection of different parameter combinations. The degree to which preferential link creation can be captured by the ODE depends on the parameter values. As expected when all link activation–deletion rates are not type dependent (i.e. *α*_{SS}=*α*_{SI}=*α*_{II} and *ω*_{SS}=*ω*_{SI}=*ω*_{II}), the agreement between ODE and simulations is excellent as illustrated in figure 7 for a range of different parameter values. In this regime, the assumption of independence of links of the triple closure performs even better than on a static network as the random breaking and creating of links helps to dissipate some of the high correlation in triplets (e.g. [*III*]) observed with *SIS*-type dynamics. Random activation–deletion helps to decrease the building of correlation in the network that may otherwise invalidate the independence assumption. Good agreement between the two models also holds for parameter values that are much different from the equal activation and deletion rates scenario. Figure 8 shows that agreement can also be found for more realistic parameter values where, for example, efforts to control a disease may require a higher *ω*_{SI} rate compared with other link deletion rates. In such situations however, care has to be taken since the agreement depends on whether the precise parameter values will lead to moderate correlations or correlations that can be captured by the ODE model. Finally, in figure 9, the bistability regime is illustrated for both the simulation and the ODE model. We note that the networks resulting from the coupled dynamics on the upper branch have on average six nodes that become isolated. As expected, the ODE cannot capture this but, we can take this into account by decreasing population size from *N*=100 to *N*=94. This results in excellent agreement between the two models. The plots in figure 9 have been obtained via continuation method where the steady state obtained at a value *τ*_{0} of the transmission rate is used as the initial condition for a new set of simulations with a smaller value of the transmission rate, i.e. *τ*_{0}−*δτ*. This approach has also been used with the full set of ODEs with time-dependent average connectivity *n*(*t*). The qualitative agreement is good and highlights that the pairwise model can be used as an exploratory tool when mapping out full system behaviour. This is especially useful where models with many parameters make this exploration almost impossible via simulations alone.

## 6. Discussion

The paper set out to carry out a systematic analysis of a model where network dynamics is coupled with a simple disease dynamics model. Network dynamics is based on link activation and deletion that can depend on link type. A step-by-step approach has been taken where the network dynamics has been modelled and studied in isolation, without disease dynamics, and then in the presence of node labelling that was stationary in time. The key result is the analysis of the pairwise ODE model and its comparison with the simulation where both network and disease dynamics act concurrently. Different modelling approaches have been proposed and discussed, and their performance relative to each other has been investigated. The models studied have ranged from exact stochastic models based on Kolmogorov equations and simple low-dimensional ODE models to simulation often with good agreement between complex simulations and ODE models. We have highlighted that approximate ODE models are desirable as they are more tractable compared with simulation and can be used as a more rigorous tool to investigate the full spectrum of behaviour, especially when faced with more complicated models with a large number of parameters.

Pairwise models have been developed in the context of static networks with applications in epidemiology and ecology (Matsuda et al. 1992; Keeling 1999; Rand 1999). This approach has been generalized to dynamic networks coupled with simple epidemics and has been shown to have the potential to be a useful modelling tool that complements simulation models and aids analysis. However, many open problems remain regarding the validity of these simple ODE models when compared with full simulation. To date, there is no coherent framework with theoretical results where, for example, the convergence, in some appropriate limit, of a stochastic process to a deterministic model is formalized. This holds even for the static network case where the need for such developments has been highlighted (House & Keeling 2011). In our particular case, there are many parameter regimes where agreement between the ODE model and simulation is good. These include, for example, regimes where either of the two dynamics is fast or slow compared with the other. This time-scale separation reflects that when considered separately, the ODE model tends to perform better. Furthermore, if the parameters are such that the system is dominated by nodes in either susceptible or infected state, then the agreement is again good given that many of the link activation and deletion processes are dormant. Hence, the dynamics of the network is predominantly governed by the creation and deletion of a unique link type, i.e. non-preferential link dynamics.

As previously highlighted, and also confirmed here, the interaction of the two dynamical process can lead to networks that are either highly clustered, display some specific structure or are fragmented in sub-networks dominated by *S* or *I* nodes even with some being completely isolated. Obviously in such cases, the current pairwise model will fail to approximate the simulation correctly but alternatives models, more sophisticated pairwise or other type of models that keep track of the number of links or even the type of neighbours a node has (Lindquist et al. 2011) may prove to give a more satisfactory result. Such models could better capture the correlations driven by preferential link dynamics, heterogeneity in node degree or other node, pair of larger scale property, that cannot be captured by the basic pairwise model.

The results of the paper also highlight the importance of coupling network and disease dynamics. This is easy to motivate from a practical viewpoint and our model illustrates that this is a worthwhile exercise with findings that complement the body of knowledge focusing on dynamics on static networks. As demonstrated in the paper, the network dynamics can be exploited as a control mechanism aimed at stopping or reducing the impact of an epidemic and there remains much work to understand how the link-type-dependent cutting and creation as well the relative time scale of the two processes impact on the type and nature of observable system behaviour. The proposed network dynamic model is basic and, while it sheds light on some interesting and relevant behaviour, it needs to be adapted based on some empirical evidence of social interaction and perturbations to it by adversities such as epidemics. While the applications of the modelling framework was in terms of epidemic models, coupling network dynamics with dynamical processes unfolding on the network remains a very current topic in areas such as computational neuroscience (notably in developmental neuroscience, e.g. Hagmann et al. (2010) and more generally the problem is recognized in Bullmore & Sporns (2009)) and we foresee that this will lead to various network dynamics and coupling mechanisms with important contributions in the future on both the modelling and analysis front.

## Acknowledgements

Istvan Z. Kiss acknowledges support from EPSRC (EP/H001085/1). Timothy Taylor is funded by a PGR studentship from MRC, and the Departments of Mathematics and Informatics at Sussex University. Péter L. Simon acknowledges support from OTKA (grant no. 81403) and from the European Union and the European Social Fund (financial support to the project under the grant agreement no. TÁMOP-4.2.1/B-09/1/KMR).

- Received May 31, 2011.
- Accepted December 21, 2011.

- This journal is © 2012 The Royal Society