Random graph models for dynamic networks
Abstract
We propose generalizations of a number of standard network models, including the classic random graph, the configuration model, and the stochastic block model, to the case of timevarying networks. We assume that the presence and absence of edges are governed by continuoustime Markov processes with rate parameters that can depend on properties of the nodes. In addition to computing equilibrium properties of these models, we demonstrate their use in data analysis and statistical inference, giving efficient algorithms for fitting them to observed network data. This allows us, for instance, to estimate the time constants of network evolution or infer community structure from temporal network data using cues embedded both in the probabilities over time that node pairs are connected by edges and in the characteristic dynamics of edge appearance and disappearance. We illustrate our methods with a selection of applications, both to computergenerated test networks and realworld examples.
I Introduction
Networked systems, such as social, technological, and biological networks, have been the subject of a vigorous research effort over the last decade Newman10 (), but most work has focused on static networks that do not change over time. In reality, almost all networks do in fact change, with nodes or edges appearing or disappearing over time, and a body of new work aimed at quantifying, modeling, and understanding such temporal or dynamic networks has recently emerged, driven in part by the increasing availability of relevant data HS12 (); Holme15 ().
Data on dynamic networks comes in a variety of forms, but the most common form, and the one we consider in this paper, is that of a set of snapshots of network structure taken at successive times, usually (though not always) evenly spaced. Such sets are a special case of a more general “multilayer” or “multiplex” network, meaning a set of different networks defined on the same set of nodes Boccaletti14 (); DGPA16 (). Multiplex networks include many nondynamic kinds, such as social networks with different types of interactions between the same set of actors. Our focus in this paper, however, is solely on dynamic networks. We also limit ourselves to networks defined on a fixed and unchanging set of nodes, so that only edges appear and disappear, not nodes. Our goal is to show how some of the most fundamental models for static networks can be generalized to the dynamic case and to demonstrate how comparisons between these models and realworld data can help us better understand the structure of the data.
Our models are built upon the assumption that the appearance and disappearance of network edges obeys a continuoustime Markov process. That is, edges appear and disappear by making transitions from present to absent or vice versa with fixed rates per unit time. Crucially, however, these rates can differ from edge to edge and, as we will see, they can have quite complex structure. If the rates of appearance and disappearance of edges are low compared to the rate at which we observe our snapshots of network structure, then consecutive snapshots will be correlated, a crucial feature of real dynamic networks. In a friendship network, for instance, one expects to still be friends next week with most of the same people one is friends with this week, so there is a strong effect of previous friendship on future friendship probability, which is reproduced by our models. In simple analyses of dynamic networks, researchers have in the past treated snapshots as independent measurements of network structure, analyzing each snapshot separately using conventional static network methods Holme15 (). This, however, ignores the often strong correlations between snapshots and thereby also ignores a potential rich source of information hidden in the data. For instance, it could be the case that two links in the network are each present in half the snapshots, but that one of these links flickers on and off rapidly, while the other one turns on and off more slowly; our Markov process model would distinguish these links as clearly different, while the more traditional model of analysis, ignoring correlations, would not.
An alternative way to think about our approach is that the fundamental unit of analysis in our calculations is not a single network but the entire history of a network, and hence that the appropriate models are those that generate entire histories. These are the models that we study in this paper.
Within this class of dynamic network models, we show how to formulate dynamic equivalents of the classic random graph, the configuration model, and the widely used stochastic block model, specifically its degreecorrected variant. As we will show, dynamic variants of simple random graph models, for instance, allow us to properly define and measure the probability of an edge or the degree of a node in an evolving network, things that otherwise present a difficult moving target, and to attach values to the rate parameters that quantify how quickly edges appear and disappear. Dynamic generalizations of block models allow us to infer largescale structure, including (but not limited to) community structure, using maximum likelihood methods akin to those developed previously for the static case.
A number of other authors have previously considered dynamic generalizations of basic network models, particularly the stochastic block model XFS10 (); Yang11 (); KL13 (); MM15 (); Ghasemian15 (); Xu14 (); MRV15 (); HXA15 (); SSTM15 (); XH13 (). The ordinary static version of the stochastic block model divides network nodes into groups or communities and then places edges between them with probabilities that depend on group membership. Dynamic variants of this idea have been investigated in which nodes can change their community membership over time, which can cause edge probabilities also to change and hence edges to appear or disappear from one snapshot to the next. Versions of this idea include the dynamic mixedmembership model of Xing et al. XFS10 () and the multigroup membership model studied by Yang et al. Yang11 () and Kim and Leskovec KL13 (). In Matias and Miele MM15 () and Ghasemian et al. Ghasemian15 (), group memberships can change but edges at successive times are independent conditioned on the groups. Xu Xu14 () has studied a dynamic block model with edge dynamics controlled by a Markov process, which has some elements in common with our approach. Matias et al. MRV15 () have considered “longitudinal” networks where contacts between nodes are governed by a Poisson process.
A little further from our focus in this paper are the multilayer stochastic block models studied for instance in Refs. HXA15 () and SSTM15 (). As with dynamic models, these models generate a set of different networks or “layers” built upon the same set of nodes, but there is now no ordering of the layers or any assumption that adjacent layers are more similar than distant ones. Han et al. HXA15 () have used such multilayer models to derive more consistent estimation of community structure for certain data sets than those derived from standard stochastic block models. Stanley et al. SSTM15 () studied a variant in which different layers (“strata” in their terminology) are generated from different underlying parameters.
In Section II we lay out the general principles behind our models, giving definitions and a variety of mathematical results for each of our models in turn. In particular, we describe dynamic versions of three static models: the Erdős–Rényi random graph, the configuration model of random graphs with a specific degree sequence, and the degreecorrected stochastic block model. We also provide efficient algorithms for statistical inference using these models, showing how to perform a maximumlikelihood fit of each one to observed data. In Section III we apply these models and algorithms to synthetic (i.e., computergenerated) test networks and to realworld examples, including technological and social networks. In Section IV we briefly describe our conclusions.
Ii Dynamic network models
Each of the models we study has a fixed number of nodes, plus edges between them that appear and disappear as the network evolves over time. Starting from some initial condition at time , our models generate continuoustime network histories, where edges appear and disappear at a sequence of realvalued times. In some data sets, events like these can be observed directly, for instance in a network of telephone calls where we are given the time and duration of each call. Here, however, we assume that the network is only observed at a set of further snapshots, evenly spaced at integer times . Including the initial state there are, thus, a total of distinct snapshots. Note, however, that the network is assumed to exist and to continue to evolve unobserved between the snapshots.
The fundamental idea behind all of the models we consider is that the edge between each node pair obeys a continuoustime Markov process, appearing and disappearing with constant rates, though the rates can differ from one node pair to another, depending on various latent properties of the nodes. By choosing this dependence appropriately, we can model various kinds of dynamic network structure, including fluctuating density, degree distribution, or community structure.
To make our discussion more concrete, consider a particular pair of nodes in the network. Let us define to be the rate (in continuous time) at which an edge appears between these two nodes where previously there was none, and let us define to be the rate at which an existing edge disappears. If we denote by and respectively the probabilities that there is and is not an edge between our nodes at time then
(1)  
(2) 
and hence satisfies the master equation
(3) 
which has the solution
(4) 
where is an integration constant and we have made use of .
Now suppose that there is no edge between our two nodes at time , i.e., that , which corresponds to the choice . Then the probability of having an edge between our nodes at the next snapshot of the network, at time , is equal to , which takes the value
(5) 
This is the probability of appearance of an edge between one snapshot and the next. Similarly we can show that the probability of disappearance of an edge is
(6) 
It will be more convenient to define our models in terms of persnapshot probabilities such as these, which can always be calculated if necessary from the fundamental rates and .
ii.1 The dynamic random graph
The random graph , famously studied by Erdős and Rényi in the 1950s and 60s ER59 (); ER60 (), is perhaps the most fundamental of all network models. In this model edges are placed between nodes pairs independently with probability (or not with probability ). In this section we define the first and simplest of our dynamic network models as a direct dynamic counterpart to the random graph.
The definition of the model is straightforward. Starting from some initial state at time , at every snapshot each node pair not connected by an edge at the previous snapshot gains an (undirected) edge with probability , or not with probability . Similarly each existing edge disappears with probability or not with probability . The net result after timesteps is a sequence of snapshots which can be represented by a set of symmetric adjacency matrices having elements if nodes and are connected by an edge in snapshot and otherwise.
In the limit of long time , the average probability of an edge between two nodes in this model is given by Eq. (4) to be , the same for every node pair. Hence the stationary distribution of the model is simply the random graph . It is in this sense that the model is a dynamic generalization of the random graph.
This is a particularly simple example of the class of models we study—we will look at more complex ones shortly—but even so there are various reasons to be interested in a model of this kind. One could use it for instance to compute the time variation of network properties such as connectivity or component sizes, or the density of specific subgraphs—computations akin to the classic calculations of Erdős and Rényi and others for the static case ER59 (); ER60 (). Our primary interest in this paper, however, is in the use of this and other models as tools for understanding observed network data, using methods of statistical inference: we fit the model to the data by the method of maximum likelihood and the parameters of the fit tell us about our data in much the same way that fitting a straight line through a set of points can tell us about the slope of those points.
Suppose that we have a set of observed snapshots of some network, measured at uniform intervals over time. If we hypothesize that the data were in fact generated from our dynamic random graph model, then the probability, or likelihood, that we observe this particular set of snapshots, given the parameters of the model, has the form
(7) 
Note that we have separate terms in this expression for the first snapshot and all succeeding snapshots. The first snapshot differs from the others because it has no preceding snapshots and hence its probability is not conditioned on those before it. The probabilities of all later snapshots, on the other hand, depend on the preexisting state of the network. Because of the assumption that network evolution follows a Markov process, each snapshot only depends directly on the immediately preceding snapshot, hence the inclusion of in the second product.
The two probabilities and are straightforward to write down. The first, which represents the probability of observing given no information about the previous history of the network, is equal to the stationary probability of an edge or nonedge within the model, which as we have said is for an edge, or for a nonedge. Hence
(8) 
The second probability is only a little more complicated, taking one of four values for edges that appear or not and ones that disappear or not:
(9) 
Substituting (8) and (9) into Eq. (7) then gives us the full likelihood for our data. In fact, as is often the case, it is more convenient to work with the logarithm of the likelihood, which has its maximum in the same place. Taking the log of (9), we have
(10) 
Given the likelihood, we can estimate the parameters and by maximizing, which gives
(11)  
(12) 
Note that these expressions differ from the naive estimates of and , given by the number of times an edge appeared or disappeared divided by the number of times it could potentially have done so. The difference arises because the initial state of the network is chosen from the stationary distribution, and the probability that in this initial state itself depends on and . As the effect of the initial state becomes progressively diluted relative to the effect of the other snapshots and Eqs. (11) and (12) converge to the naive values.
Because appears on the righthand side of (11) and (12), calculating the rates and requires us to find selfconsistent solutions to the equations. In fact, it is possible to eliminate the dependence on on the righthand side and derive explicit closedform equations, but the expressions are somewhat complicated. In practice we have found it simpler just to solve Eqs. (11) and (12) by iteration from a suitable initial condition.
What do these equations tell us? For a given data set, they give us an optimal estimate—better than the naive estimate—of the rate at which edges appear and disappear in our network. This gives us information about the correlation between adjacent snapshots. The combined values of and also give us the maximumlikelihood estimate of the average density of the network, via the average probability of an edge.
This model, however, while illustrative, is not, in practice, very useful. Like the static random graph which inspired it, it is too simple to capture most of the interesting structure in real networks, and in particular it generates networks with Poisson degree distributions, wholly unlike those of realworld networks, which typically have broad and strongly nonPoisson distributions. In the world of static network models, this latter shortcoming is remedied by the configuration model, a more sophisticated random graph that can accommodate arbitrary degree distributions MR95 (); NSW01 (). In the next section, we show how to define a dynamic equivalent of the configuration model.
ii.2 Dynamic random graphs with arbitrary expected degrees
The configuration model is a model of a random graph with a given degree sequence MR95 (); NSW01 (). One fixes the degree of each node and then places edges at random subject to the constraints imposed by the degrees. This can be achieved in practice by endowing each node with “halfedges” and choosing a matching of halfedges uniformly at random from the set of all possible matchings. In the limit the expected number of edges falling between nodes and in this model is , where is the total number of edges in the network, and the actual number of edges between each pair of nodes is Poisson distributed with this mean. There is nothing in this model to stop a pair of nodes having two or more edges connecting them—a socalled multiedge—and in general there will be some multiedges in networks generated using the configuration model. Selfloops—edges connecting a node to itself—can and do also appear. Although this is not realistic behavior for most realworld networks, versions of the configuration model that explicitly forbid multiedges and selfloops are much harder to work with than those that do not. Moreover, if the degree distribution has finite mean and variance, the expected number of multiedges and selfloops in the network is constant, independent of , so they have vanishing density as . For these reasons, one normally puts up with the presence of a few multiedges and selfloops for the sake of simplicity.
A commonly studied variant of the configuration model, which is easier to treat in some ways, involves explicitly placing between each node pair a Poissondistributed number of edges with mean . In this variant, sometimes called the Chung–Lu model after two of the first authors to study it CL02b (), the numbers of edges between node pairs are independent random variables, making analysis simpler. The price one pays for this simplicity is that the degrees of individual nodes are no longer fixed, themselves being Poissondistributed (and asymptotically independent) with mean . Thus in this case represents not the actual degree but the expected degree of a node. (The random graph of Erdős and Rényi, with mean degree , is then the special case of this model where for all .)
In this section we define a dynamic analog of the Chung–Lu model in the sense of the current paper: its edges have a dynamics chosen so that the stationary distribution of the model is precisely the Chung–Lu model. Since the Chung–Lu model can contain multiedges, we consider a process for adding and removing edges slightly different from the one of the previous section, such that each pair of nodes can have any nonnegative number of edges connecting it. Specifically, for each node pair we consider the Poisson process where edges are added at rate , and each of the existing edges is removed independently at rate . Thus is incremented with rate , and decremented with rate .
Let denote the probability that a node pair has edges at time . Then satisfies the master equation
(13) 
We can solve this equation by defining a generating function , multiplying both sides of (13) by , and summing over to get
(14) 
The general solution to this equation is
(15) 
where is any oncedifferentiable function of its argument satisfying , the latter condition being necessary to fulfill the normalization requirement for all .
In the limit of long time we have , which is the generating function of a Poissondistributed variable with mean . Hence the number of edges between any pair of nodes in this model is Poissondistributed in the stationary state. If we make the choice
(16) 
for some set of values , with as previously and any value of , then the mean number of edges between nodes and is . In other words, the stationary state of this model is precisely the Chung–Lu model with expected degrees .
This then defines our model: to generate a dynamic network with nodes, we specify the expected degree for each node and the parameter . We generate the initial state of the network from the Chung–Lu model with these expected degrees, and then generate future states by adding edges between each node pair at rate and removing existing edges at rate . We sample snapshots of the resulting network at integer intervals which, along with the initial state at , comprise the total snapshots generated by the model. We represent these snapshots by adjacency matrices .
One could use this model for various purposes, such as making calculations of expected structural properties, but our principal interest here is again in fitting the model to observed network data. As before we achieve this by maximizing a likelihood function, which has the same basic form as previously:
(17) 
The first probability on the righthand side is straightforward to write down, since we know that the stationary distribution places a Poissondistributed number of edges between nodes and with mean . Thus
(18) 
which is independent of .
The second probability is more involved, but the calculation is simplified by noting that even though the model can possess multiedges, the observed network data will normally have at most a single edge between any pair of nodes, so that the only allowed edge transitions are the appearance and disappearance of single edges.
Suppose that a given node pair is connected by zero edges at time . Then, setting in Eq. (15), we find that , which implies that one timestep later at we have
(19) 
where we have made use of Eq. (16) and for convenience defined the quantity
(20) 
which (by analogy with our use of the same symbol in Section II.1) is equal to the total probability that an existing edge disappears during a single unit of time, i.e., between two successive snapshots.
The probabilities and of a transition from zero edges to, respectively, zero or one edges in a single timestep are then equal to the probabilities and of having zero or one edges at . These are given by the zeroth and first coefficients in the expansion of in powers of :
(21)  
(22) 
By a similar method we also have
(23)  
(24) 
where we have ignored terms of second and higher order in in (24).
We can now write down the transition probability as a function of :
(25) 
Substituting this into Eq. (17) and taking logs, we get the following expression for the loglikelihood in our model:
(26) 
where
is the total number of newly appearing edges in the observed data, and similarly
(27) 
Then, differentiating (26) with respect to , we find that the optimal value of is the positive solution of the quadratic equation
(28) 
Similarly, differentiating with respect to and bearing in mind that , we find that obeys
(29) 
which has the solution
(30) 
The sum in this expression is the number of edges initially connected to node plus the number that later appear. The divisor is the effective number of independent measurements of an edge that we make during our snapshots. If , so that edges never appear or disappear, then in effect we only have one measurement of each edge—the initial snapshot at . Conversely, if , so that every observed edge immediately disappears on the next snapshot, then all snapshots are independent and the number of independent measurements is . Thus Eq. (30) measures the number of observed edges between node pairs divided by the number of independent observations of each node pair.
Equations (28) and (30) give us the maximumlikelihood estimates the rate parameter and the expected degrees of the nodes. We note two points:

These equations have to be solved selfconsistently, since the first equation depends on via and the second depends on .

Neither nor are equal to their naive estimates from the data. One might imagine, for instance, that would be given by the average of over all snapshots, but our results indicate that the maximumlikelihood estimate differs from this value.
Both of these effects arise, as in the previous section, because of the information provided by the initial state. Because the initial state is drawn from the stationary distribution, which depends on the model parameters, we can make a better estimate of those parameters by taking it into account than not. On the other hand, the advantage of doing so dwindles as becomes large and vanishes in the limit.
We could use these equations, for example, to define in a principled fashion an equivalent of the “degree” for a node in a dynamic network. The actual degree of a node in such a network is a fluctuating quantity, but using our results one can define a single number for each node that, like the degree in a static network, is a measure of the propensity of that node to connect to others. We give some examples in Section III.
ii.3 Dynamic block models
The stochastic block model is a random graph model of a network that incorporates modules or community structure—groups of nodes with varying densities of within and betweengroup edges. The standard stochastic block model, first proposed by Holland et al. in 1983 HLL83 (), is the communitystructured equivalent of the random graph of Erdős and Rényi, but like the latter it has shortcomings as a model of realworld networks because the networks it generates have Poisson degree distributions. The degreecorrected stochastic block model KN11a () is a variant on the same idea that is analogous to the model of Chung and Lu CL02b (), allowing us to choose any set of values for the expected degrees of nodes, while also generating a communitystructured network. In this section we define a dynamic equivalent of the degreecorrected block model along similar lines to the models of previous sections and show how it can be used to infer community structure from dynamic network data.
The standard (static) degreecorrected block model divides a network of nodes into nonoverlapping groups labeled by integers . Let us denote by the group to which node belongs. Then we place a Poissondistributed number of edges between each node pair with mean equal to , where is a degreelike parameter and is a further set of parameters which control the density of edges within and between each pair of groups. If the diagonal elements are greater than the offdiagonal ones, this model generates networks with conventional “assortative” community structure—dense ingroup connections and sparser betweengroup ones—although other choices of are also possible and are observed in realworld situations.
This description does not completely fix the parameters of the model: they are arbitrary to within a multiplicative constant, since one can multiply all the in any group by a constant and divide the same constant out of without affecting the behavior of the model. This is why we refer to as a “degreelike parameter”—it plays a role similar to degree in the configuration model, but is arbitrary to within a groupdependent multiplicative constant. Following KN11a (), we remove this ambiguity by making a specific choice of normalization, that the sum of within any group should be 1:
(31) 
where is the Kronecker delta. This gives us constraints, one for each of the groups, and hence fixes all the remaining degrees of freedom.
To generalize this model to the dynamic case we again divide our nodes into groups and assign to each of them a degreelike parameter satisfying (31). We generate an initial state drawn from the static degreecorrected block model with these parameters. We then generate a history for the network by adding edges between each node pair at rate
(32) 
and removing existing edges independently at rate , where and are respectively the groups to which and belong. Note the similarity between Eqs. (16) and (32), the primary differences being that the parameter now depends on the group memberships and that the factor has been replaced by the quantity , which also depends on the group memberships. By the same argument as before, the number of edges between and in the stationary state is Poisson distributed with mean
(33) 
which makes the stationary state of this model equivalent to the degreecorrected stochastic block model as desired.
Also by the same argument as before, we can calculate the transition rates for edges to appear and disappear between one snapshot and the next, which are
(34)  
(35)  
(36)  
(37) 
Here
(38) 
is the total probability for an existing edge between nodes in groups and to disappear in the unit of time between successive snapshots. (Also as before we have in Eq. (37) discarded terms beyond leading order in the small quantities .)
By fitting this model to observed network data, we can determine the parameters , , and , along with the group assignment parameters . The likelihood as a function function of the four sets of parameters , , , and takes the form
(39) 
The first probability on the right is straightforward, taking the value
(40) 
by definition (which is independent of ), while the second can be expressed in terms of the transition probabilities, Eqs. (34) to (37). The resulting expression for the loglikelihood is
(41) 
where
(42) 
and
(43) 
which is the total number of edges that appear between groups and in the observed data. Similarly,
(44)  
(45) 
Differentiating Eq. (41) with respect to now gives us
(46) 
and differentiating with respect to gives a quadratic equation again:
(47) 
(Note that in order to perform the derivatives correctly, one must take into account the fact that and , although it turns out that the end result is the same as would be derived by naive differentiation, ignoring these equalities.)
Differentiating (41) with respect to , and normalizing appropriately, gives us
(48) 
The selfconsistent solution of Eqs. (46), (47), and (48), now gives us the parameters of the model.
If we want to convert the degreelike parameter into a true degree, we can do this by noting that the expected degree of node in the stationary state of this model is equal to the sum of the expected number of edges between and every other node, which is
(49) 
where we made use of Eq. (31) in the final equality. Hence the degrees are simply proportional to , with a constant of proportionality that can be easily calculated once we have the values of from Eq. (46).
This still leaves us to calculate the maximumlikelihood estimates of the group assignments . To do this, we substitute our estimates of the parameters back into the loglikelihood, Eq. (41), to get the socalled profile likelihood, which is then maximized over the group assignments . Note that there is no need to calculate the last term in the likelihood since, by Eq. (46), it is equal to , which is independent of the group assignments and hence has no effect on the position of the maximum.
Maximization of the profile likelihood over the values of is harder than maximizing with respect to the other parameters, since the values of the are discrete. We perform the maximization numerically, using a heuristic algorithm analogous to that used for the static block model in KN11a (), which was in turn inspired by the classic Kernighan–Lin algorithm for graph partitioning KL70 (). Starting from a random group assignment, we move a single node to a different group, choosing from among all possible such moves the one that most increases (or least decreases) the profile likelihood. We repeat this process, making a chain of successive singlenode moves, but with the important qualification that each node is moved only once. When all nodes have been moved once, we reexamine every state passed through during the process to find the one with the highest profile likelihood, then take that state as the starting point for a new repetition of the same algorithm. We continue repeating until no further improvement in the profile likelihood is found. As with many other optimization algorithms, the results can vary from one run to another because of the random initial condition, so one commonly performs several complete runs with different initial conditions, taking as the final answer the output of the run that gives the highest overall value of the profile likelihood.
An alternative way to fit our model would be to use an expectation–maximization (EM) algorithm in which the model parameters are assigned their maximumlikelihood values but one computes an entire posterior distribution over divisions of the network into groups. The latter distribution, being a large object, is normally evaluated only approximately, either by Monte Carlo sampling or using a belief propagation algorithm DKMZ11a () in which nodes pass each other estimates of their (marginal) probabilities of belonging to each group. A belief propagation algorithm was used previously for a different dynamic block model in Ghasemian15 (), where each node sends messages both along “spatial” edges to its neighbors in each snapshot and along “temporal” edges to its past and future selves in adjacent snapshots. A similar approach could work in the present case, although our model differs from that of Ghasemian15 () in assuming unchanging group memberships but correlated edges where Ghasemian15 () makes the opposite assumption of timevarying group memberships but independent edges between snapshots.
Iii Applications
In this section we give examples of fits of dynamic network data to the dynamic configuration model of Section II.2 and the dynamic block model of Section II.3.
iii.1 Synthetic networks
Our first set of examples make use of synthetic data sets—computergenerated networks with known structure that we attempt to recover using the maximumlikelihood fit. We demonstrate this approach using the dynamic block model of Section II.3 and the test networks we use are themselves generated using the same model. We look in particular at the case where the expected degree parameters for all nodes are the same, equal to a constant . For the tests reported here we use . At the same time we varying the strength of the community structure, encapsulated in the parameters , according to
(50) 
Here is diagonal (all elements with are zero), is a flat matrix (all elements are the same), and is an interpolating parameter. Thus by varying we span the range from a uniform random graph with no community structure () to a network in which all edges lie within communities and none between communities (), so that the communities are completely disconnected components.
We similarly vary the rate constants according to a second parameter , also lying in , such that
(51) 
which interpolates between values that are the same for all groups and the heterogeneous choice , which can be anything we choose. Note that while varying does not change the expected degree or average density of edges in the network, it does change how rapidly edges appear and disappear. Thus controls the extent to which the dynamics of the network, as opposed to merely its average behavior, gives additional information about the community structure.
Once the parameters are fixed, we generate a set of networks, which in our tests have nodes divided into two groups of equal size. For each network we generate an initial state followed by up to five further snapshots. The initial state is generated from the stationary distribution (i.e., from a traditional degreecorrected block model) and the following snapshots are generated according to the prescription of Section II.3.
We now apply the fitting method of Section II.3 to these networks to test whether it is able to successfully recover the community structure planted in them. Success, or lack of it, is quantified using the normalized mutual information DDDA05 (); Meila07 (), an informationtheoretic metric that measures the agreement between two sets of group assignments. As traditionally defined, a normalized mutual information of 1 indicates exact recovery of the planted groups while 0 indicates complete failure—zero correlation between recovered and planted values.
Figure 1 shows the results of our tests. In panel (a) we fix , so that is uniform and block structure is indicated only by the relative abundance of edges within and between groups. We use a value of , meaning that 40% of extant edges disappear at each timestep. The different curves in the figure show the normalized mutual information as a function of the parameter which measures the strength of the community structure, for different numbers of snapshots from to . As we can see, our ability to recover the planted structure diminishes, and eventually fails completely, as the structure becomes weaker, but this effect is partly offset (as we might expect) by increasing the number of snapshots—the more snapshots we use the better we are able to infer the community structure. For larger numbers of snapshots, the algorithm is able to surpass the known “detectability threshold” below which community detection is impossible for single, static networks DKMZ11a (), which is indicated by the vertical dashed line in the figure. In other words the algorithm is able to integrate information about the network over time in order to better determine the shape of the communities.
In Fig. 1b we set , so that , choosing the value of to be along the diagonal and off the diagonal, meaning that withingroup edges are somewhat more persistent—more likely to be conserved from one snapshot to the next—than betweengroup edges. This behavior provides another signal of community structure, in addition to the differing timeaveraged edge probabilities, which the algorithm can in principle use to determine group memberships. And indeed the results of Fig. 1b reflect this, showing that the algorithm is able to determine group memberships even well below the detectability threshold, but only when is large. If is small, then it becomes difficult to determine the values of from the data, and hence difficult to determine group membership for small . This point is discussed further below.
In Fig. 1c we fix and vary between zero and one using values as previously, and , . With there is now no signal whatsoever of community structure present in the positions of the edges. The only clue to the group assignments lies in the rate of appearance and disappearance of edges within and between groups. As we would expect, the algorithm is unable to identify the communities at all when or , but as grows for the algorithm assigns a larger and larger fraction of nodes to the correct groups, with better performance for larger values of . These results suggest the existence of a new detectability threshold as a function of , with location tending to zero as . (A threshold like this was observed, for instance, by Ghasemian et al. Ghasemian15 () in their model, discussed in Section II.3, which has a transition as a function of both the strength of community structure and the relevant rate parameters.)
iii.2 Realworld examples
We have also tested our models against a number of empirical data sets representing the structure of realworld dynamic networks. We give three examples representing networks drawn from technological and social domains, finding in each case that our dynamic models and their associated algorithms perform better than static methods.
Internet graph:
Our first example is a network representation of the structure of the Internet at the level of autonomous systems (ASes), the fundamental units of global packet routing used by the Internet’s Border Gateway Protocol. The structure of the Internet changes constantly and is well documented: a number of ongoing projects collect snapshots of the structure at regular intervals and make them available for research. Here we use data from the CAIDA AS Relationships Dataset CAIDA (), focusing on four snapshots of the network’s structure taken at threemonth intervals during 2015. The spacing of the snapshots is chosen with an eye to the rate of growth of the network. The Internet has grown steadily in size over the several decades of its existence, and it is still growing today, but this growth is not captured by our models. To ensure better fits, therefore, we first restrict our data to the set of nodes that are present in all of our snapshots, and second choose snapshots that span a relatively short total time. Thus our four snapshots were chosen to be sufficiently far apart in time that the network sees significant change between one snapshot and the next, but close enough that the size of the network does not change significantly.
We fit our Internet data to the dynamic version of the configuration model described in Section II.2, which gives us a way to determine the parameter that controls the rate of appearance and disappearance of edges as well as the effective degrees of nodes in the network. For the rate parameter we find a maximum likelihood value of , which indicates a fairly slow rate of turnover of the edges in the network. Recall that is the average probability that an edge will vanish from one snapshot to the next, so this value of implies that over 90% of edges remain intact between snapshots. As discussed in Section II.2, one could make a naive estimate of the rate at which edges vanish simply by counting the number that do, but that estimate would be less accurate than the maximumlikelihood one.
Our fit also gives us estimates of the degree parameters from Eq. (30). Figure 2 shows a histogram of the frequency distribution of estimated degrees for the Internet derived in this manner. Again, we could make naive estimates of the degrees, for instance by assuming snapshots to be independent and averaging the raw degrees of their nodes across snapshots. This would be a correct estimator of the in the limit of a large number of snapshots, meaning it will tend to the correct answer eventually, but it would be less than ideal. In particular, our estimate of the error on the values it gives would be wrong. By assuming the snapshots to be independent, we effectively assume that we have more measurements than we really do and hence underestimate the variance. For instance, if we observe that the naive degree of a node is unchanging for many snapshots in a row, we may conclude that the average of those values has a very small statistical error, because the fluctuations are small. This, however, would be erroneous if the small fluctuations are actually just a result of the fact that the network is only changing rather slowly.
Error estimates are not the only thing that will be affected by improperly using a naive degree estimate. The values of the degrees themselves can also be affected if the snapshots are strongly correlated, which they are in this case because of the small value of . Strongly correlated snapshots will tend to give a node the same or similar degree on successive snapshots, but Eq. (30) implies that in this case our estimate of should actually decrease over time (as becomes larger in the denominator while the numerator remains constant). A naive estimate on the other hand would remain unchanged. At first sight the decrease in the maximumlikelihood estimate may appear counterintuitive, but it has a simple physical interpretation: for a node that truly has a constant value of , we would expect additional edges to appear occasionally, at a rate dependent on the value. If we do not see any edges appearing, therefore, it implies our initial estimate of the degree was too high and we should revise it downward.
The maximumlikelihood estimator can, on the other hand, also have problems of its own if the the model we are fitting is not a perfect description of the data. In the case of the Internet we see two possible sources of disagreement between data and model. First, even though the number of nodes in the network is held fixed, the number of edges is observed to grow over time—the network is becoming more dense. This effect is not included in our model, which assumes constant expected density. Second, we see some evidence that the removal of edges is not uniform as our model assumes, but that edges connected to highdegree nodes disappear at a higher rate than those connected to lowdegree ones. Both of these behaviors could potentially affect our results. (It is interesting to ask whether and how the model could be generalized to include them, though we leave pursuit of these questions for future work.)
Friendship network:
Our second example focuses on a set of social networks from a study by Michell and West of friendship patterns and behaviors among school students in the UK MW96 (). Highschool students at a school in the west of Scotland were polled about their friendship patterns, each student being allowed to name up to twelve friends, and they were also asked about their drinking, smoking, and drug use habits. The entire exercise was conducted a total of three times, at yearly intervals, with the same group of students. The study looked at all students in the school, but the most detailed data were collected for a subset of 50 girls within the larger population and it is on this subset that we focus here.
The researchers were interested in the extent to which substance use behaviors correlated with friendship patterns. They found that although there was no single factor that would completely explain the friendships of the students, the network of friendships did display homophily according to substance use, meaning that students with similar use patterns were more likely to be friends PW03 (); PSS06 ().
In our analysis, we divide the students into three groups: those who do not drink, smoke, or take drugs on a regular basis; those who exhibit one of these three behaviors; and those who exhibit two or more. We then ask whether it is possible to detect this division into groups based on network structure alone, without any knowledge of student behaviors. We find that when using the dynamic version of the degreecorrected block model described in Section II.3 it is indeed possible to determine the groups, and to do so with better accuracy than can be achieved by standard static methods. Specifically, we compare results from our dynamic block model to those from the static degreecorrected block model fitted to an aggregate network formed from the union of the three snapshots.
Figure 3 shows three pictures of the overall aggregate network of friendships. Each picture is laid out identically, but with different coloring. In panel (a) the three colors represent the ground truth, with green, yellow, and red denoting students who engaged in zero, one, or two or more of the behaviors studied respectively. Panel (b) shows the communities found in the network by fitting to the dynamic block model. Though not perfect, this fit places 64% of the nodes in their correct groups. A random coloring, for comparison, would get only 33% right. Panel (c) shows the results from the standard static algorithm applied to the aggregated network. This fit places only 52% of the nodes in their correct groups.
Proximity network:
Our third example is another social network, a network of physical proximity between students in a high school in France MFB15 (). The data were collected using electronic proximity detectors worn by the participants. The detectors record the presence and identity of other detectors in their vicinity at intervals of 20 seconds. The data were collected over five consecutive days, but only a half day’s worth of data were collected on the last day, which we discard, leaving four full days to work with. We construct one snapshot for each day and consider there to be an edge between two participants in a snapshot if three or more contacts between them are recorded during the relevant day. Requiring a minimum number of contacts in this way helps to remove spurious signals from the data, as discussed in SBPT16 (). We also restrict our study to those nodes that are present in all snapshots.
The students in the study were divided among three subject specialties: mathematics/physics, physics/chemistry, and engineering. Each specialty was further divided into three classes, so there are a total of nine classes in the data. We attempt to recover these classes from the network data alone, without other information, using both the dynamic model of this paper and a traditional static degreecorrected block model applied to the aggregated network. In this case both methods do well, which is perhaps unsurprising, given that the edges within each group are significantly denser than those between groups. Figure 4 shows the results for the dynamic model in panel (a) and the static model in panel (b). As we can see, both models achieve good classification of the nodes into their classes, though the dynamic model performs slightly better. The error rate—the fraction of incorrectly labeled nodes—is 4.1% for the dynamic model of panel (a) and 5.7% for the static model of panel (b).
The primary benefit of the dynamic model in this case, however, lies not in its ability to recover the communities but in what it reveals about the dynamics of the network. In addition to the communities themselves, the dynamic model also returns values for the rate parameters that can reveal features of the data not seen in the simple static fit to the aggregate network. Of particular interest in this case are the parameters , which measure the relative rates at which edges change within and between groups. Our fit gives estimates of
(52) 
In other words, connections are not only more likely between participants in the same class or specialty, but they are also more persistent, in some cases by a wide margin—only about 6% of connections persist from one snapshot to the next between individuals in the different specialties for example, but almost 50% persist within classes. (There is some variation in values of among classes and specialties; the results above are only an approximate guide based on average values for each type.)
Iv Conclusions
In this paper we have introduced dynamic generalizations of some of the bestknown static network models, including the Erdős–Rényi random graph, the configuration model, and the degreecorrected stochastic block model. We have also derived and implemented efficient algorithms for fitting these models to network data that allow us to infer maximumlikelihood estimates of rates of change, node degrees, and community structure. We have tested the performance of our models and algorithms on synthetic benchmark networks as well as on a selection of data sets representing realworld examples of dynamic networks.
There are a number of directions in which this work could be extended. First, we have focused exclusively on edge dynamics here, but there are also networks in which nodes appear and disappear and it would be a natural generalization to study the dynamics of nodes also, or of both edges and nodes together. (We could also allow node properties, such as expected degrees or community memberships, to change over time, as some other authors have done.) Second, the assumption of continuoustime Markov processes for the edge dynamics is a particularly simple one, which could be relaxed to encompass more complicated situations. Third, in our community detection calculations we assume we know the number of communities the network contains, but in many cases we do not have this information. Methods have been developed for determining the number of communities in static networks and it is an interesting question whether those methods can be extended to the dynamic case as well.
Acknowledgements.
The authors thank Aaron Clauset for useful conversations. This research was supported in part by the US National Science Foundation under grants DMS1107796 and DMS1407207 (MEJN) and by the Army Research Office under grant W911NF12R0012 and the John Templeton Foundation (CM).References
 M. E. J. Newman, Networks: An Introduction. Oxford University Press, Oxford (2010).
 P. Holme and J. Saramäki, Temporal networks. Phys. Rep. 519, 97–125 (2012).
 P. Holme, Modern temporal network theory: a colloquium. The European Physical Journal B 88(9), 1–30 (2015).
 S. Boccaletti, G. Bianconi, R. Criado, C. I. del Genio, J. GomezGardenes, M. Romance, I. SendinaNadal, Z. Wang, and M. Zanin, The structure and dynamics of multilayer networks. Physics Reports 544, 1–122 (2014).
 M. De Domenico, C. Granell, M. A. Porter, and A. Arenas, The physics of multilayer networks. Preprint arxiv:1604.02021 (2016).
 E. P. Xing, W. Fu, L. Song, et al., A statespace mixed membership blockmodel for dynamic network tomography. The Annals of Applied Statistics 4, 535–566 (2010).
 T. Yang, Y. Chi, S. Zhu, Y. Gong, and R. Jin, Detecting communities and their evolutions in dynamic social networksâa bayesian approach. Machine Learning 82, 157–189 (2011).
 M. Kim and J. Leskovec, Nonparametric multigroup membership model for dynamic networks. In Advances in Neural Information Processing Systems, pp. 1385–1393 (2013).
 C. Matias and V. Miele, Statistical clustering of temporal networks through a dynamic stochastic block model. Preprint arXiv:1506.07464 (2015).
 A. Ghasemian, P. Zhang, A. Clauset, C. Moore, and L. Peel, Detectability thresholds and optimal algorithms for community structure in dynamic networks. Phys. Rev. X 6, 031005 (2016).
 K. S. Xu, Stochastic block transition models for dynamic networks. Preprint arXiv:1411.5404 (2014).
 C. Matias, T. Rebafka, and F. Villers, Estimation and clustering in a semiparametric poisson process stochastic block model for longitudinal networks. Preprint arXiv:1512.07075 (2015).
 Q. Han, K. Xu, and E. Airoldi, Consistent estimation of dynamic and multilayer block models. In Proceedings of The 32nd International Conference on Machine Learning, pp. 1511–1520 (2015).
 N. Stanley, S. Shai, D. Taylor, and P. J. Mucha, Clustering network layers with the strata multilayer stochastic block model. Preprint arXiv:1507.01826 (2015).
 K. S. Xu and A. O. Hero III, Dynamic stochastic blockmodels: Statistical models for timeevolving networks. In Social Computing, BehavioralCultural Modeling and Prediction, pp. 201–210, Springer (2013).
 P. Erdős and A. Rényi, On random graphs. Publicationes Mathematicae 6, 290–297 (1959).
 P. Erdős and A. Rényi, On the evolution of random graphs. Publications of the Mathematical Institute of the Hungarian Academy of Sciences 5, 17–61 (1960).
 M. Molloy and B. Reed, A critical point for random graphs with a given degree sequence. Random Structures and Algorithms 6, 161–179 (1995).
 M. E. J. Newman, S. H. Strogatz, and D. J. Watts, Random graphs with arbitrary degree distributions and their applications. Phys. Rev. E 64, 026118 (2001).
 F. Chung and L. Lu, The average distances in random graphs with given expected degrees. Proc. Natl. Acad. Sci. USA 99, 15879–15882 (2002).
 P. W. Holland, K. B. Laskey, and S. Leinhardt, Stochastic blockmodels: Some first steps. Social Networks 5, 109–137 (1983).
 B. Karrer and M. E. J. Newman, Stochastic blockmodels and community structure in networks. Phys. Rev. E 83, 016107 (2011).
 A. Decelle, F. Krzakala, C. Moore, and L. Zdeborová, Inference and phase transitions in the detection of modules in sparse networks. Phys. Rev. Lett. 107, 065701 (2011).
 B. W. Kernighan and S. Lin, An efficient heuristic procedure for partitioning graphs. Bell System Technical Journal 49, 291–307 (1970).
 L. Danon, J. Duch, A. DiazGuilera, and A. Arenas, Comparing community structure identification. J. Stat. Mech. 2005, P09008 (2005).
 M. Meilă, Comparing clusterings—an information based distance. Journal of Multivariate Analysis 98, 873–895 (2007).
 The CAIDA AS relationships dataset. http://www.caida.org/data/asrelationships.
 L. Michell and P. West, Peer pressure to smoke: The meaning depends on the method. Health Education Research 11, 39–49 (1996).
 M. Pearson and P. West, Drifting smoke rings. Connections 25, 59–76 (2003).
 M. Pearson, C. Sieglich, and T. Snijders, Homophily and assimilation among sportactive adolescent substance users. Connections 27, 47–63 (2006).
 R. Mastrandrea, J. Fournet, and A. Barrat, Contact patterns in a high school: A comparison between data collected using wearable sensors, contact diaries and friendship surveys. PLOS One 10, e0136497 (2015).
 J. C. Silva, L. Bennett, L. G. Papageorgiou, and S. Tsoka, A mathematical programming approach for sequential clustering of dynamic networks. The European Physical Journal B 89, 1–10 (2016).