Draft updated August 26, 2016

Formal analysis of the emergent structural properties of dynamic networks is largely uncharted territory. We focus here on the properties of forward reachable sets (FRS) as a function of the underlying degree distribution and edge duration. FRS are defined as the set of nodes that can be reached from an initial seed via a path of temporally ordered edges; a natural extension of connected component measures to dynamic networks. Working in a stochastic framework, we derive closed-form expressions for the mean and variance of the exponential growth rate of the FRS for temporal networks with both edge and node dynamics. For networks with node dynamics, we calculate thresholds for the growth of the FRS. The effects of finite population size are explored via simulation and approximation. We examine how these properties vary by edge duration and different cross-sectional degree distributions that characterize a range of scientifically interesting normative outcomes (Poisson and Bernoulli). The size of the forward reachable set gives an upper bound for the epidemic size in disease transmission network models, relating this work to epidemic modeling (Ferguson 2000, Eames 2004).

Forward Reachable Sets]Forward Reachable Sets: Analytically derived properties of connected components for dynamic networks B. Armbruster, L. Wang, M. Morris] Benjamin Armbruster
Northwestern University
and  Li Wang, Martina Morris
University of Washington

1 Introduction

There is growing interest in diffusion over dynamic networks (a recent review may be found in [\citenameHolme & Saramäki, 2012]), but formal analysis of the emergent structural properties of dynamic networks is largely uncharted territory. In this paper we focus on deriving properties of forward reachable sets (FRS), a natural temporal extension of the static concept of a connected component. The FRS is defined as the set of all nodes that can be reached from an origin node via a forward reachable path over some period of time; a forward reachable path is a sequence of temporally ordered edges that connects two nodes within a specified time interval [\citenameMoody, 2002].

Our work is motivated by applications in epidemiology, where an infectious disease passes from one person to another in a population through a sequence of partnerships that form, have duration, and dissolve. Partnerships are different than “contacts” in this context; they represent repeated contact with the same person over time, and this temporal clustering of contacts can affect transmission dynamics. The duration of partnerships influences the stability of the underlying transmission network structure, and the cross-sectional degree distribution (the number of partners at one instant) influences how connectivity emerges over time. Long-term monogamous partnerships can isolate dyads, reducing epidemic potential, but if partnerships turn over rapidly or people have more than one partner at a time, epidemic potential increases. The impact of these pair formation features on epidemic outcomes has been recognized in the sexually transmitted infection (STI) literature for decades [\citenameDietz & Hadeler, 1988, \citenameMorris & Kretzschmar, 1997, \citenameKretzschmar & Dietz, 1998, \citenameBauch & Rand, 2000, \citenameFerguson & Garnett, 2000, \citenameLeung et al., 2014].

The FRS can be thought of as the upper bound for the transmission process on a partnership network: an idealized infectious disease for which the probability of transmission is 1, and transmission is instantaneous. We model the FRS as an emergent structural property of the dynamic network, generated by micro-level processes that govern node and link dynamics: the mean cross-sectional degree , the edge dissolution rate , and the node exit rate . Conditional on these parameters, we vary the cross-sectional degree distribution, comparing Bernoulli to Poisson. These two distributions are chosen because they map to the normative continuum that guides partnership dynamics. The Bernoulli, at one end, represents serial monogamy, where only one partner at a time is allowed so each partnership is entirely dependent on the presence of another, while the Poisson at the other end assumes all partnerships are formed and dissolved independently, so there is no restriction on the number of partners one can have at the same time. Populations with the same mean number of partners per person can vary along the continuum defined by these two distributions. Populations with a highly skewed partnership distribution, such as the power-law distribution, are not represented in our analysis. However, long-tailed distributions are usually found for cumulative degree (e.g., number of partners in the last year) [\citenameHamilton et al., 2008], rather than cross-sectional degree. We will discuss the implications of long-tailed distributions in Section 3.1 and in the discussions (Section 7.3).

Our analysis focuses on the effect of the degree distribution on the properties of the FRS, conditional on mean degree, edge duration, and node exit rates. Working in a stochastic framework, we find analytic solutions for the threshold for growth of the FRS and the mean and variance of the active FRS size over time. Our analysis considers three phases in the evolution of the FRS: the initial exponential growth phase where we can make the approximation of an infinite size network, the logistic growth phase where the finite network effect kicks in, and the equilibrium where the active FRS is either extinct or varies around a non-zero size. We evaluate the stochasticity in the FRS by calculating the probability of extinction, and the variance of the size of the FRS in the initial growth phase.

Focusing on the FRS allows us to greatly simplify the pair-formation and moment closure models found in the epidemics literature because we do not track the neighbors’ status in our equations. As a result, we are able to obtain closed-form analytic solutions to several traditional epidemic potential indicators, as a function of basic parameters that can be obtained from egocentrically sampled networks in standard sample surveys. The methods we develop here can therefore be used to empirically investigate epidemic potential with easily collected network data.

2 Definitions and Setup

2.1 Precise definitions and assumptions

A dynamic network can be represented as a list of edges with activities over intervals (Holme 2012): , where each 4-tuple contains the two endpoints of the edge (), the edge formation time , and the dissolution time . The active intervals () of an edge must have positive length (), and may not overlap another interval on the same edge. The networks in our analysis are undirected.

We also model node entry and exit (analogous to birth and death in a dynamic population model). Once a node exits the network, all of its edges dissolve and it can no longer form edges; the exit is permanent. Each node has only one active lifetime, but each edge may be associated with multiple active intervals. The two endpoints of an active edge must both be active nodes.

The network is assumed to be homogeneous and in equilibrium. All nodes have the same rate of exit , regardless of their degree. All nodes with degree have the same edge formation rate . All active edges have the same rate of dissolution . All rates are constant across time. All dynamics are assumed to follow a memoryless stochastic process with the specified rates, which determine how quickly the network is changing. The memoryless process implies that partnership durations are exponentially distributed.

Note that an edge dissolution may be caused by the death of one of the two endpoints, or by a separation process at rate . We can write the dissolution rate as and the separation rate as . Since the separation rate must be positive, there is a constraint that .

We examine dynamic networks with two different cross-sectional degree distributions: Poisson(k) and Bernoulli(p). The Poisson degree distribution corresponds to a uniform edge formation process, where for any node degree, resulting in an Erdős-Rényi random network at equilibrium when the number of nodes is infinite. This type of graph is also referred to as a “simple” or Bernoulli random graph in different literatures. The Bernoulli degree distribution restricts edge formation to nodes that have degree 0, meaning if and for . It corresponds to a strictly monogamous social network constraint. Let be the mean degree of a cross-sectional network taken at equilibrium. Equivalently, can be the mean instantaneous degree of a node at equilibrium.

Figure 1: Forward reachable path: the black horizontal lines represent nodes, and the blue shaded regions represent active edges. The red line is a forward reachable path from node 1 to 4, within . There is a path even though the cross sectional network at has no active connecting edges. The dotted path also connects node 1 to 4, but takes a longer time.

Given the set of node and edge activities, we can define a forward reachable path as a time-ordered sequence of active edges between an initial and a destination node, within a time interval. An equivalent definition using cross-sectional networks can be found in Nicosia et al.\shortcitenicosia_components_2012. Our definition emphasizes the endpoints of the activity intervals, which are the events that form the basis of our analysis.

  • There exists a missingforward reachable path from node to itself in any time interval (reflexivity). There exists a forward reachable path from to in the time interval if there exists a node such that the following two conditions hold:

    • There is a path from to in the time interval , with .

    • There is an active edge where and .

In our definition for the forward reachable path, there is no delay between lateral moves in the path, going from one node to another. In applications like disease transmission or airline connections, an incubation or transit period can be specified by requiring that .

Unlike connectivity on a static network, the forward reachable path is always directed, even if the network is undirected. Connectivity via forward reachable path is not symmetric: existence of a path from node to does not imply a path from to . However, it is transitive if we account for the time coordinate: if there is a path from to in the time interval and a path from to in the time interval , then there is a path from to in the interval .

  • The missingforward reachable set (FRS) given is the set of all nodes where a forward reachable path exists between and in the time interval . The active forward reachable set is the subset of nodes in the FRS that are still active at .

One immediate consequence of the reflexivity in forward reachable path is that if there exists a path between and in the interval , then there exists a path between those two nodes for all intervals with . So, the forward reachable set (FRS) is non decreasing in time. However, the active FRS can get smaller when nodes exit the network, and a node that has exited the network will not be able to form any more edges, so it will not participate in the growth of the FRS. When all nodes in the active FRS have exited, the active FRS becomes extinct, and its size will stay at zero.

The dynamics of the network (formation, dissolution, and death) is not influenced by whether a node is in the FRS, but the FRS is completely determined by the dynamic network. So, the FRS can be viewed as a structural property of the dynamic network, analogous to component size in a static network.

2.2 Goals for analysis

Our goal is to model the size of the forward reachable set over time on a random dynamic network, as a function of the network parameters: mean cross-sectional degree , cross-sectional degree distribution (Bernoulli vs Poisson), edge dissolution rate , and node exit rate (Table 1). We want to define appropriate approximations and find analytic solutions for the mean and variance of the FRS size, and compare the results to simulation results.

The networks in the model are finite size, undirected, with population dynamics (node entry and exit). The number of active nodes and edges of the network are assumed to be in equilibrium so that there is a balance between the node entry and exit rates, and between the edge formation rate and the edge dissolution rate . This rate balancing will be the foundation of our analysis (Section 3).

The outcome variables of interest are the expected size and variance of the active part of the FRS , as a function of time . If there is no node exit , then is equivalent to the FRS size. The variation is considered over all initial nodes in the network at . There are three phases in the evolution of the FRS: the initial exponential growth phase where we can make the approximation of an infinite size network, the logistic growth phase where the finite network effects kick in, and the equilibrium where the active FRS either becomes extinct or varies in size around a non-zero value.

We split our analysis into three sections. For the initial exponential growth, we make the approximation that a randomly selected node is already in the FRS with probability 0; this is equivalent to assuming an infinite network size. Given this assumption we can derive exact solutions for the expectation and variance of the FRS size (section 4). For the logistic growth phase, we modify our equations to include the finite network size effect, using a mean-field approximation (section 5). For the equilibrium persistence, we use the mean-field approximation to derive the expected fraction of the network in the active FRS () (section 6). We use simulations to show the probability of extinction in this stochastic process, and derive this probability analytically. Our main results will compare the behavior of the active FRS between Bernoulli and Poisson degree distributions, for various values of mean degree, dissolution rate, and node exit rate. This comparison has strong implications for the role of concurrent partnerships in the spread of a sexually transmitted disease in a sexual network (section 7).

3 Preliminary Analysis

 Parameter Definition
Mean instantaneous degree
Edge dissolution rate, per edge
Node exit rate, per active node
Edge separation rate, per edge (derived)
Edge formation rate, per eligible node (derived)
FRS growth rate (derived)
Table 1: Parameters and notation

In this section we discuss the preliminaries for the derivation of the equations for the FRS growth in Section 4. These include calculating the size of the connected components; relating the rate of partnership formation to the network parameters; and deriving how the cumulative or “aggregate” degree of a node grows over time. A more rigorous derivation of these properties can be found in [\citenameLeung et al., 2012].

Let be a random variable denoting the size of the connected component attached to a randomly chosen node, on a cross-sectional (static) network. For a Bernoulli degree distribution, clearly . For Poisson \shortcitenewman_component_2007, when


We do not consider the case since it would lead to a giant connected component, with . These results are valid on a homogeneous network of infinite size. For finite or non-homogeneous networks, the component size can be obtained via simulation [\citenameMorris & Kretzschmar, 2000].

The total dissolution rate is a combination of partnership separation and death of either partner, so , and the expected partnership duration is . The separation rate is .

Next we derive the partnership formation rates per eligible node, and for the Poisson and Bernoulli cases respectively, given the average degree , edge dissolution rate , and node exit rate . We focus on the degree over time of a random active node, and model as a Markov birth-death process on in the Poisson case and on in the Bernoulli case. Since the degree is only relevant when the node is alive, we are implicitly conditioning on that, so the death rate of the birth-death process (not the death rate of the node) is : an active node may lose a partner due to separation at rate , or partner death at rate .

In steady-state, , balancing the birth and death rates of the birth-death process. For the Poisson case where all nodes are eligible for formation, this balance is . For the Bernoulli case, only the degree-0 nodes are eligible to form an edge, and only degree-1 nodes can dissolve an edge. Noting that fraction of nodes with degree 1 is , and the fraction of nodes with degree 0 is , the equilibrium balance is and therefore, .

3.1 Cumulative degree

As a final piece of the preliminaries we look at the aggregate or cumulative degree. Often, surveys do not ask “how many partners do you have at this moment?” but commonly “how many partners (or new partners) have you had in [time period]?”. The first is , the equilibrium distribution of . Here we examine on the second. Let be the cumulative number of partners of some node from its birth until a duration later. In the Poisson case, the node acquires partners independent of its current number of partners, so is just a Poisson process with rate . The Bernoulli case is slightly trickier since the time between new partners (i.e., when increases) is . Note that the expectation of this time is still as in the Poisson case, but the variance of this time is now which is smaller than the variance in the Poisson case, . Essentially, is a renewal process, specifically an equilibrium renewal process if . Thus,


in both the Poisson and Bernoulli cases and


in the Poisson case and asymptotically in the Bernoulli case,


using standard renewal results. Note that both grow linearly with time, the expectations are the same for Poisson and Bernoulli, and the variance is smaller for the Bernoulli case.

For the Poisson case, the number of partners within a time interval would have a Poisson distribution, if the node is actively acquiring partners for the whole interval. If however the node became active during this period, then the number of partners within the time interval would be less. This heterogeneity would lead to a more skewed distribution for cumulative degree in our model.

One related calculation that is simple to perform is the lifetime number of partners, . In the Poisson case, the number has a Geometric distribution:


due to standard properties of these distributions. In either the Poisson or Bernoulli case we can use iterated expectations to show that . Thus,


assuming births have an equilibrium number of partners. Poisson and Bernoulli degree networks have the same expected lifetime number of partners.

3.2 Simulation setup

We will compare our analytic results to simulations on random dynamic networks, both to verify the analysis and to show differences between the theory and the actual stochastic process. The networks we simulated have 500 nodes over a range of parameters , , and . The random dynamic networks are simulated as continuous time Markov processes with three types of events, using the stochastic simulation algorithm [\citenameBanks et al., 2011]. A formation event will create an active edge between the sender node and a random node in the set of eligible receiver nodes (all other nodes for Poisson, all nodes with degree zero for Bernoulli). Edge separation is defined as the dissolution of an edge that’s not due to the exit of one of the nodes. A separation event deactivates a random active edge in the network. A node exit event deactivates the node and all the edges incident on it. After the network simulation is completed, we calculate the active forward reachable set () for every active node at . We used the R package “tsna” for this calculation [\citenameBender-deMoll & Morris, 2015].

In our simulations, a node that exits the network is immediately replaced by a new node (“reincarnation”). The new node enters the network with degree 0. This is significant for a social network interpretation, but it adds an inhomogeneity to the dynamic network. The new nodes do not have the same degree distribution as the rest of the network, which can affect the equilibrium size of the forward reachable set. However, new nodes do not affect the exponential growth phase: even if they enter the network with non-zero degree, it is unlikely that they are in contact with the small initial FRS at birth. In Appendix B (Section B), we will modify our model to account for the fact that the newly born nodes have degree 0.

To match the simulated network to the theoretical parameters, we adjust the rate of edge formation of each active node and the separation rate . To match the mean degree , we set the formation rate for Poisson degree distribution to , and for Bernoulli to for “eligible” nodes of degree 0. The total dissolution rate is split up into two processes, separation and node exit. The separation rate is set as .

Figure 2: Growth of the active forward reachable set (FRS) for network parameters , , , and network size 500. The thin blue lines show the FRS trajectories starting at different initial nodes on a single simulated network with Poisson degree distribution; the lighter orange lines, for Bernoulli. The circles along the bottom represent FRS extinctions, and the thick lines are for the corresponding mean FRS sizes.

Figure 2 shows realizations of the growth of the FRS over time, for different initial nodes on a single dynamic network simulation. We can separate the growth dynamics into three phases. When the FRS size is small compared to the network, it exhibits exponential growth. When the active FRS size is the same order of magnitude as the network size, the growth curve becomes logistic and eventually hits a plateau. Although there is much variation in the time for the FRS size to take-off, eventually all the realizations starting at different nodes converge to an active FRS of the same size.

4 Exponential growth phase

For the exponential growth phase, the size of the forward reachable set is small compared to the size of the network. In the epidemiology context, this means the fraction of susceptibles in the population is close to 1 [\citenameBritton & Trapman, 2014]. We will make the approximation that when an edge forms on a node in the FRS, it connects to a node not in the FRS with probability 1. This approximation is equivalent to having an infinite size network.

4.1 General ODE Model

Our general approach models the forward reachable set as a continuous time Markov chain (CTMC). Consider a vector-valued Markov jump process with state variables and transition rate matrix . We assume all vectors are column vectors by default. The Kolmogorov backward equations [\citenameKarlin, 2014] state that for any function ,


where we define


The last equality is due to the property of transition matrices.

The processes of interest to us have special structure. They only have a few different types of jumps (i.e., transitions) and the jumps are of the following form. In state we take jumps of type at a rate linear in the state, , where is a constant. Furthermore, a jump of type transitions the state from to where may be random but is required to be independent of .

With these assumptions and choosing the identity for , we can write the backward equations as


where matrix


In the next sections we specify how our problem can be written in this form deriving the equations describing how the expected size of the FRS grows over time, first for the Poisson case, then for the Bernoulli case, before finally comparing the results. Derivations and results for the dynamics of the variance of the size of the FRS are in appendix A.

4.2 Poisson degree networks

For the Poisson case, the state variables we need are , the size of the active FRS, and , the number of inactive (i.e., dead) nodes in the FRS, so in the vector notation from above, . There are two types of transitions. The first is the formation of an edge where one endpoint is in the FRS (formation). This happens at a rate of and adds a connected component to the active FRS ( ). Thus using the vector notation, and . The second type of transition is node exit (death). It happens at a rate of , and removes a single node from the FRS ( and ). Using the vector notation, and .

Thus, equation (12) gives us the following differential equations,


At time 0, and equals the size of the connected component of the initial node of the FRS. If the initial node is chosen at random, then . Node birth events are not counted here because new nodes have degree 0 and do not add to the active FRS.

These differential equations are linear with constant coefficients, and thus have simple exponential solutions:


where the exponential growth rate , with . This makes sense from a dimensional analysis perspective: since is a growth rate, it must be the product of a rate parameter such as and a function of the dimensionless constants, and .

4.3 Bernoulli degree networks

For the Bernoulli case, in addition to and , we must add a variable to keep track of the number of degree-0 nodes in the FRS. Only nodes in are eligible to form active edges. The vector of state variables is now . Our Markov process has the following types of transitions:

  • Formation. and at rate .
    , .

  • Dissolution. at rate . Here counts the number of partnerships in the FRS.
    , .

  • Death of a degree-0 node in the FRS. , , and at rate .
    , .

  • Death of a degree-1 node in the FRS. , , and at rate .
    , .

This leads to the following set of differential equations:


As in the Poisson case, the equation for decouples from the rest; none of the other equations rely on . Thus we end up with a two dimensional linear system with constant coefficients, where the growth rate of the solution is determined by the largest eigenvalue of the coefficient matrix. The eigenvalues are . Since and , the growth rate is:

where again . Thus and the same for . In fact in both the Poisson and Bernoulli cases, the differential equation for is the same and thus, as . This growth rate has simple edge cases: for and for .

4.4 Key results

Since and are both dimensionless, dimensional analysis tells us that we can write the growth rate in the form . Equivalently, we can arbitrarily choose and plot as a function of for various choices of . Figure 3 shows this. From this we see that the Poisson growth rates go to infinity as , because most nodes are connected as part of a giant component in the cross-sectional network. The Bernoulli degree network has a max component size of 2, so its FRS growth is limited by the edge dissolution rate. We can show that the growth rate for Poisson degree distributed network is faster than Bernoulli, , given any parameter values.

Figure 3: Growth rate, . Black lines are for the Poisson case, blue lines for the Bernoulli case. Here and we show . For other values, scale and accordingly.
Lemma 4.1

for all parameter values.


Let . Then where . Now , proving the claim.

In simulations, the exponential growth rate of the expected FRS size matches up well with the analytic results Figure 2. The growth rate is higher for larger cross-sectional mean degree and larger values of the dimensionless parameter , which is proportional to the number of formations per node lifetime (Section 3). We can make the comparison between Poisson and Bernoulli degree distributions for different values of and . Figure 4 shows that the FRS growth rate is always larger for a Poisson degree network compared to a Bernoulli degree network, given the same parameter values. A Poisson degree network can achieve the same growth rate with a smaller mean degree, or fewer formations per node lifetime.

Figure 4: Comparison of the FRS growth rate across simulated networks with different parameter values. The size of the circles represent the mean FRS size at , which we can us as a proxy measure of the exponential growth rate. We only show circles for FRS sizes 3%. The two lines highlight the theoretic contours of equal growth rate.

We can define a growth threshold, , that separates the parameter values where the FRS grows exponentially from those where it shrinks or doesn’t grow. From dimensional analysis we can write the threshold as for some function or equivalently to write it in a form similar to the standard form for . In the Poisson case, , so is equivalent to . For the Bernoulli case we have . The thresholds are similar for small because for small values of , the left hand side of the Bernoulli threshold is , matching the Poisson threshold. In addition Lemma 4.1 of course implies that the Bernoulli threshold is always greater than the Poisson one (i.e., implies ). Figure 5 plots the two thresholds as functions of and the dimensionless quantity .

Figure 5: Growth threshold. Above the threshold, growth is positive. There is a range of parameters (the blue shaded region) with positive growth for Poisson, but not for Bernoulli degree networks. Note is not possible (Section 2.1).

5 Logistic growth phase

When the size of the forward reachable set grows to a significant fraction (about 10%) of network size, the assumptions we used for the exponential growth phase are no longer valid. When a new edge is formed, the probability that the node in the FRS connects to a node not already in the FRS is no longer approximately 1; it is reduced to . This saturation effect is captured by the mean-field models that modify the differential equations from Section 4.

For Poisson degree networks, the differential equation becomes


For the Bernoulli case, we get


We can solve the equations above numerically to approximate the expected FRS growth in the logistic phase. The growth rate decreases as the FRS gets larger, but it never drops below 0. Figure 6 shows how the average FRS size in simulations matches up with the mean-field solutions for a 500 node network, up until the equilibrium.

Figure 6: Long-term behavior of the active forward reachable set for a simulated network with parameters , , and network size 500 (same as in Figure 2). The non-logged y-axis shows the logistic growth curve; the mean-field logistic approximations match up with the persistent FRS size. The different trajectories show the FRS growth for different starting nodes. At equilibrium, the FRS size distribution is bimodal, with the lower mode showing the extinction probability (right panel). All the trajectories that persist eventually converge to an active FRS of the same size.

Setting the right hand side to zero allows us to solve for an approximation to the equilibrium expected size. Figure 7 shows the dependence of the mean-field equilibrium on the parameters.

Figure 7: Mean-field equilibrium forward reachable set size for . Black lines are for the Poisson case and blue lines for the Bernoulli case. Lines are shown for .

6 Equilibrium phase

At equilibrium, the FRS size distribution becomes bimodal (Figure 6). For some of the initial nodes, the FRS takes off and reaches a persistent equilibrium value. For other initial nodes, the FRS shortly becomes extinct, and will remain at zero, even when the network parameters are above the threshold for growth. Our mean-field model is deterministic and does not account for extinctions of the FRS, which happens stochastically. The mean-field solutions match up with the average persistent FRS size (for the set of initial nodes that do not have their FRS go extinct).

Figure 8: The persistent FRS size, plotted as circle size, across different network parameter values at equilibrium. The band of blue circles shows the region where extinction probability is between 25% and 75%. In the lower left corner, the FRS goes extinct for most initial nodes. In the upper right, the FRS will grow to the persistent size most of the time.

6.1 Persistent FRS size

The logistic mean-field equations from Section 5 allow us to calculate the equilibrium FRS fraction for the Poisson case, by setting the growth rate to zero:


where . Note that the threshold for a positive steady state, , coincides with the threshold for a positive growth rate in the exponential growth phase, .

For the Bernoulli case, setting the differential equations (22)–(23) to 0 gives the equilibrium FRS size:


Again, the threshold for a positive persistent FRS size at equilibrium coincides with the threshold for positive FRS growth.

Lemma 6.1

iff .


. Rearranging , . So iff . Expanding, this is equivalent to proving the claim.

The persistent FRS size at equilibrium, like the growth rate, is higher for larger cross-sectional mean degree and larger values of the dimensionless parameter (the circle sizes in Figure 8).

If we compare the equilibrium behavior for Poisson degree networks versus Bernoulli in Figure 7, there is a region of parameter values between the respective thresholds where a persistent FRS is possible for Poisson degree networks, but not for Bernoulli. For a specific value of the mean degree , this region is between


which coincides with the region in Figure 5 for positive growth rates.

For all parameter values, the equilibrium persistent size should be larger for Poisson degree distributed networks. Similar to Lemma 4.1 proving that , we can prove for the mean-field models that . This is also illustrated in Figure 7 which shows the steady-state fractions, and , as functions of and . However, we do not see this in simulations, because the simulated networks are not actually homogeneous.

Lemma 6.2



Note that and . Since , , proving the claim.

The above result assumes that nodes entering the network have the same degree distribution as the rest of the network. In our simulations, nodes enter the network with degree-0, which causes the network to have a non-homogeneous degree distribution. The degree-0 entry has important interpretations for social networks, and is common for network models with population dynamics [\citenameLeung et al., 2012]. We show that there is a significant effect on the equilibrium size of the FRS. The equilibrium FRS size increases for Bernoulli networks, and decreases for Poisson networks, switching their relative positions from our homogeneous model (Figure 6 and 9). We make adjustments to our model and derive the FRS outcome in Appendix B.

Figure 9: Comparing the persistent FRS size from the homogeneous model (Section 5) vs the adjusted degree-0 entry model (Appendix B). For network parameters , , , and network size 500 (same simulations as in Figure 2), the homogeneous model shows higher equilibrium FRS size for Poisson degree networks than Bernoulli. However, in the degree-0 entry model, and in our simulations (thin lines), the equilibrium prevalence for both are nearly the same.

6.2 Stochastic behavior and extinction

We can model the FRS extinction as a branching process, a technique common in epidemiology. Let the population be the nodes in the active FRS. Let the “offspring” of a node be the non-FRS nodes it contacts, during the time interval when it is in the active FRS. Because stochastic extinctions mostly occur within the first few generations, the FRS is small compared to the network, as in the exponential growth phase (Section 4). So, we can approximate the number of offspring by the cumulative degree of a node, . This approach is similar to the concept of basic reproduction number in epidemiology [\citenameLeung & Kretzschmar, 2015].

Let be the probability generating function (PGF) of the offspring distribution . Then the probability of extinction is the smallest non-negative root of the equation [\citenameKarlin, 2014]


There are two parts to a node’s cumulative degree : the initial neighbors at the time it enters the FRS, and the number of formations over its lifetime . The initial degree has a distribution specified by the model (either Poisson or Bernoulli). From Section 3.1 we know that for Poisson degree distributions, the cumulative number of formations has a geometric distribution (7).

So, the PGF of is


where . We can solve this equation numerically to find the extinction probability for different values of the parameters for a Poisson degree network. As expected, the extinction probability goes to 1 for parameter values below the growth threshold (Section 4.4). We lack the distribution of lifetime number of formations for Bernoulli degree networks for this calculation.

For parameter values below the threshold, the FRS will go extinct for most initial nodes. For high values of and , the FRS will go to persistence most of the time. But in a band of parameter values just above the threshold, the FRS exhibits high variability, where there is large uncertainty in the equilibrium outcome (Figure 8).

Poisson Bernoulli
Edge separation rate
Edge formation rate
Expected lifetime partners
Initial FRS growth rate
Threshold for growth
Equilibrium prevalence
Extinction probability , solve for
Table 2: Closed-form expressions for our key results. The expressions depend only on the network parameters (mean degree), (edge dissolution rate), and (node exit rate). For simplicity, we defined and . The equilibrium prevalence is shown for the base model, where new nodes are indistinguishable from existing nodes. The prevalence in our simulations, where new nodes enter with degree 0, is derived in Appendix B.

7 Discussion

We developed analytic expressions and simulation tools for exploring the properties of forward reachable sets in stochastic dynamic networks with open populations. Our model provides insight into the thresholds, rates of growth, and equilibrium states of the FRS, representing these as emergent features of more basic network properties: degree means and distributions, edge dissolution rates, and node exit rates. All of these are properties that can be measured in sample surveys, which allows our methods to be grounded in empirical data when desired. We find that the two degree distributions compared here, Bernoulli and Poisson, produce quantitative and qualitative differences in the properties of the FRS when conditioning on all other parameters. The epidemic threshold is lower and the growth rate is faster in Poisson networks than in Bernoulli networks, for a wide range of underlying mean degrees and partner acquisition rates. There is also a broad band of parameter values around the threshold for both network types where the variability in the growth rate and extinction process is maximal. In this band, the possibility of extinction leads to a bimodal equilibrium FRS size distribution, one that emerges at the beginning of the process, not just over the long term.

A natural application of these findings is in epidemiology, where the FRS represents the maximal epidemic potential in a population. There is a substantial literature on the population dynamics of infectious disease, and a smaller but growing literature on the impact of partnership network structure on these dynamics. Most of this work relies on deterministic models and mean field approximations for analytic results. Our work therefore contributes to this literature in several ways. First, our stochastic analysis describes the variance of the size of the FRS (maximal disease prevalence) on a non-trivial network over time. Second, we derive analytic solutions for the expected growth rate and growth threshold for this stochastic process. Third, we are able to derive a mean field approximation for the equilibrium prevalence from our stochastic process of pair formation and dissolution with both vital dynamics and non-trivial network structure; and compare this to stochastic simulations. And fourth, under these conditions we derive an analytic expression for the extinction probability for the Poisson degree distribution case. The methods we present are specifically designed to capture the impact of stochastic variability in a systematic way, so that its practical implications can be evaluated.

There are many different types of infectious diseases; HIV, in particular, motivated our work, and this influenced several aspects of the analysis. We note some of the resulting implications and limitations below.

7.1 Disease States

Since nodes do not leave the FRS and nodes cease forming ties after they exit, this model is analogous to an “SI” (susceptible-infected) epidemic, and the findings are not relevant for diseases with recovery, either back into susceptible status (SIS) or with immunity (SIR). Our model is appropriate for HIV, but not for other diseases such as measles or chlamydia.

7.2 Parameter Values

We examined a range of parameter values that are both plausible for sexual partnership networks and bound the transitional region for the FRS. The mean active lifetime is set at 40 years; 35-40 is often used in HIV simulation analyses [\citenameBaggaley et al., 2006, \citenameHallett et al., 2008, \citenameEaton & Hallett, 2014]. Mean relationship duration ranges from 2.5 to 15 years, and the exponential distribution results in both a large fraction of shorter relationships as well as a right-skewed exponential tail [\citenameLeung et al., 2012]. Mean cross-sectional degree ranges from 0.2 to 0.6 partners, similar though a bit lower than that reported in a study of young adults in the U.S. [\citenameMorris et al., 2009]. Combining these parameters, the average lifetime number of partners ranges from 0.5 to 10, similar to the range reported in many surveys, both in the U.S. and elsewhere [\citenameMorris, 1993, \citenameMorris et al., 2010, \citenameHamilton & Morris, 2010].

7.3 Degree Distribution Comparison

In the context of sexual partnership networks, the Bernoulli degree distribution represents a restrictive rule of serial monogamy for partnerships, while the Poisson instead allows individuals to have multiple concurrent partners. While these distributions do not represent any specific population, they bracket a reasonable range for general heterosexual populations. Our results provide additional support to the growing literature that finds concurrent partnerships increase both epidemic potential and variability [\citenameGoodreau et al., 2012, \citenameMorris & Kretzschmar, 1997, \citenameEaton et al., 2011, \citenameLeung & Kretzschmar, 2015].

There is extensive literature documenting a highly skewed, power-law-like cumulative degree distribution in sexual partnership networks, reviewed in [\citenameHamilton et al., 2008]. However, the degree distributions in our model refer to the number of partners at one moment in time, rather than cumulative over the last year. Although our model would not reproduce a power-law network, we can get closer to a long-tailed distribution by relaxing the homogeneity assumptions in Section 2.1. Specifically if we allow each node to form connections at different rates (but still constant over time), the network will have more nodes with a large cumulative degree, given the same network density. The impact of such heterogeneity on the FRS will be a topic of further investigation.

7.4 Additional Heterogeneity

A realistic model would need to account for a wide range of additional demographic details that influence partnership network structure and dynamics: group-specific mean-degree, dissolution, and mortality rates, and differences in behavior based on disease status. For example, if nodes in the FRS have a reduced partnership formation rate, or an increased mortality rate, then the FRS growth rate will obviously decrease, and the threshold for growth would shift toward larger parameter values for mean degree and dissolution rate. If these demographic details only affect the dynamic network, our model can account for it by adding more compartments, at the cost of additional equations [\citenameFerguson & Garnett, 2000], and the loss of closed-form results.

Additionally, a memoryless stochastic process (“Markov assumption”) does not account for different types of relationships. In real networks, partnership durations are often bi-modal (long and short term relationships), and partnership durations and type may both depend on the presence of other partnerships for each node. As such, our framework seems better suited for gaining broad intuition into the impact of network structure and dynamics on reachability, than for making detailed predictions for specific epidemics.

7.5 Partnership status at entry

In the finite population models and in our simulations, we assume that new nodes enter the network as non-FRS with degree zero (Section 3.2). As a result, nodes in the FRS have higher mean degree than nodes outside. In the context of disease spread, this would mean that infected individuals will on average have more contacts than the non-infected. The implications are derived in Appendix B. For Bernoulli networks (i.e., under the assumption of serial monogamy), the implication is that the fraction of singles inside the FRS is lower, especially when the FRS size is large. As a result, a node in the FRS has a higher chance of pairing with an outside node, which leads to a higher rate of discordant contact, and shifts the equilibrium prevalence higher. For Poisson networks, the implications is that the mean degree is lower outside the FRS, leading to smaller components outside the FRS, which reduces the FRS growth rate and shifts the equilibrium prevalence lower. This produces the outcome seen in Figure 9: a higher expected equilibrium FRS size for the Bernoulli networks. If we instead let the new nodes have the same mean degree as the existing ones, the rankings are reversed, and Poisson networks have the higher equilibrium FRS size (as we derive analytically in Lemma 6.2). We were a bit surprised that this seemingly innocuous assumption would have such an impact on a key outcome of interest. It suggests that both analytic and simulation studies of epidemic dynamics should pay attention to this assumption in order to avoid inappropriate artifacts.

7.6 Reachability vs. Epidemic Spread

Reachability implicitly defines “transmission” as incorporation into the FRS — it is instantaneous with probability 1 on first contact between discordant nodes (i.e., an infinite transmission rate). This represents the theoretical maximum any transmission process could reach on a specific network. A model for the spread of any specific pathogen would require finite, possibly heterogeneous and possibly time-varying transmission rates across discordant pairs. For example, a model for HIV transmission would need to capture variation in transmission probability by stage of infection, demographic subgroup, treatment status, and the use of prophylaxis.

Some of the implications of a finite transmission rate are obvious: it would increase the probability of epidemic extinction and the parameter values needed to cross the persistence threshold, and would reduce the growth rate and equilibrium prevalence. There is no reason to expect that the qualitative results for the Bernoulli vs. Poisson comparison would be changed — a Poisson network will always offer multiple pathways for transmission, and larger connected components for any fixed parameter value set — but the range of parameter values in which this difference matters will change to reflect the new region of epidemic persistence.

Deriving analytic results for stochastic models with a finite transmission rate, however, will be more complicated. The epidemic growth rate is proportional to the number of discordant edges, and one would now need to keep track of these over time until transmission actually occurs. Deterministic models that keep track of discordant edges, like the pairwise moment closure model [\citenameEames & Keeling, 2002], pair-formation model [\citenameKretzschmar & Dietz, 1998], and edge-based compartmental model [\citenameMiller et al., 2011], offer analytic approaches and some solutions for the growth rate, threshold, and final prevalence of an epidemic. However, these models often assume that the network is static, that the degree distribution is binomial (equivalent to our Poisson model for large networks) or regular (everyone has the same number of partners), that triad closure is negligible, and all of them ignore stochasticity. Triad closure is clearly important for finite transmission rates since any node can only be infected once in an SI process, and we have already discussed the importance of the other factors.

7.7 Stochasticity of Epidemic Outcome

Most previous papers on transmission processes did not provide any insight into the stochastic variability of epidemic outcomes, or did so for a restrictive model (for example, [\citenameDangerfield et al., 2009] calculated the variance of the outcome on a static regular network). When infections (like HIV) have a slow transmission rate or networks (like sexual partnership networks) have slow link turnover, then the stochasticity of the initial epidemic growth will be important, and the epidemic outcome cannot be cleanly partitioned by the growth threshold. For parameters above the threshold, there is still a chance that an epidemic can go extinct. Conversely, when the parameters are below the threshold, there is a chance that an epidemic can persist for a long time.

Our findings suggest that this stochastic variability is important, especially in the band around the persistence/growth threshold. The practical implications are that one may find very different dynamics and outcomes in populations with nearly identical network structure and dynamics. This variability makes it more difficult to infer epidemic potential from observed epidemic outcomes alone, and more difficult to predict epidemic outcomes going forward. Quantifying this uncertainty, in the context of a finite transmission rate, is an important subject for further investigation.

Appendix A Variance analytics

We now describe how to determine ODEs for variances, . Note that and thus, . More generally, for a vector-valued process,


From section 4 we already have equations for . We then use the approach in section 4.1 to derive equations for . Recall, for jumps of type , at rate . Thus,


This ((12) and (31)) is still a set of linear differential equations. Specifically, component of (31) is a linear combination of an extended set of state variables ,


a.1 Poisson case

Applying this framework to the transitions in the Poisson case (see section 4.2), we derive the following equations:


where as before we ignore since it does not influence the other variables. Together with the original equation for , this is a linear system of differential equations with eigenvalues and . This reasonably suggests that the standard deviation of grows at the same rate as the expectation. Specifically, we have the closed-form solutions


If we add in the equations relating to ,


then the largest eigenvalue of the system for remains as before.

a.2 Bernoulli case

We now apply the framework to the transitions for the Bernoulli case (see section 4.3 for a list of the transitions). We again ignore as it does not affect the other variables, and derive the following equations:


The largest eigenvalue of a 3x3 system involves the solution to a cubic equation and is thus a mess. Nevertheless, we conjecture that the largest eigenvalue here is less than based on numerical experiments.

For completeness, here are the equations involving (the equation for is omitted as it is the same as for the Poisson case):


Appendix B Adjustments to the Mean field Approximation

In our simulations, nodes enter the network with degree-0, which causes the network to have a non-homogeneous degree distribution. Leung et al. \shortciteleung_dynamic_2012 showed that the effect on the network is negligible if the partnership dissolution rate is fast compared to the lifetime of each individual node. However we show that there is still a significant effect on the equilibrium size of the FRS. The equilibrium FRS size increases for Bernoulli networks, and decreases for Poisson networks, switching their relative positions from our homogeneous model (Figure 9). We make adjustments to our model and derive the FRS outcome below.

b.1 Bernoulli degree networks

In a homogeneous Bernoulli degree network, the probability for a degree-0 node (“single”) in the FRS to connect with a single not in the FRS is , where is the active FRS size, and is the network size (Equation (22)). If nodes enter the network as a single (and non-FRS), the percent of singles is higher for the non-FRS group, compared to the FRS group. So, a single in the FRS has a higher chance of finding an eligible partner not in the FRS than in a homogeneous network, and we get a larger equilibrium prevalence than expected. This effect is negligible when the FRS size is small compared to the network size, so it will not affect the exponential growth phase.

We need to keep track of the number of singles not in the FRS (), and use the conditional probability in the “Single” column in Table 3, instead of the marginal probability.

 In FRS  Not in FRS
Table 3: Two-way table of relationship status vs FRS status

We add to the state variables of our Markov process, along with the new transition types in Table 4. Note that all deaths in our simulations are automatically replaced with a new birth as , to keep the network size constant.

Transition Rate Change in state (A,W,V)
WV formation (infection)
WW formation
VV formation
WW separation
VV separation
W death
V death
FRS paired death
Non-FRS paired death
Table 4: New transition types for Bernoulli degree networks

We also need an additional differential equation to track changes in . The modified equations are:

and the initial conditions are , , . Replace the state variables on the right hand side with their expectations to obtain the deterministic mean-field equations. These equations can be solved numerically.

b.2 Poisson degree networks

In a Poisson degree network, nodes connect to any other node in the network with equal probability, so we do not have the same situation as above. However, the mean degree in a network with degree-0 births is different for nodes in the FRS versus those not in the FRS. Since nodes enter the network as degree-0 and non-FRS, the mean degree for non-FRS nodes will be lower than for the overall network. For a Poisson distributed network, lower mean degree leads to lower mean component size, so we need to adjust Equation (14) to the following,

The logistic mean-field equation (21) can be adjusted in the same way. We can obtain the mean component size for non-FRS nodes, , via simulation, and solve the differential equation as before. Since this adjustment lowers the number of nodes added to the FRS in a formation event, it will lower the equilibrium prevalence for a Poisson degree network.


  • [\citenameBaggaley et al., 2006] Baggaley, Rebecca F, Garnett, Geoff P, & Ferguson, Neil M. (2006). Modelling the Impact of Antiretroviral Use in Resource-Poor Settings. Plos med, 3(4), e124.
  • [\citenameBanks et al., 2011] Banks, H. T., Broido, Anna, Canter, I, Gayvert, Kaitlyn, Hu, Shuhua, Joyner, Michele, & Link, Kathryn. (2011). Simulation algorithms for continuous time markov chain models.
  • [\citenameBauch & Rand, 2000] Bauch, C., & Rand, D. A. (2000). A moment closure model for sexually transmitted disease transmission through a concurrent partnership network. Proceedings of the royal society b: Biological sciences, 267(1456), 2019–2027.
  • [\citenameBender-deMoll & Morris, 2015] Bender-deMoll, Skye, & Morris, Martina. (2015). tsna: Tools for Temporal Social Network Analysis.
  • [\citenameBritton & Trapman, 2014] Britton, Tom, & Trapman, Pieter. (2014). Stochastic Epidemics in Growing Populations. Bulletin of mathematical biology, 76(5), 985–996.
  • [\citenameDangerfield et al., 2009] Dangerfield, C. E., Ross, J. V., & Keeling, M. J. (2009). Integrating stochasticity and network structure into an epidemic model. Journal of the royal society interface, 6(38), 761–774.
  • [\citenameDietz & Hadeler, 1988] Dietz, Klaus, & Hadeler, K. P. (1988). Epidemiological models for sexually transmitted diseases. Journal of mathematical biology, 26(1), 1–25.
  • [\citenameEames & Keeling, 2002] Eames, Ken TD, & Keeling, Matt J. (2002). Modeling dynamic and network heterogeneities in the spread of sexually transmitted diseases. Proceedings of the national academy of sciences, 99(20), 13330–13335.
  • [\citenameEaton & Hallett, 2014] Eaton, Jeffrey W., & Hallett, Timothy B. (2014). Why the proportion of transmission during early-stage HIV infection does not predict the long-term impact of treatment on HIV incidence. Proceedings of the national academy of sciences, 111(45), 16202–16207.
  • [\citenameEaton et al., 2011] Eaton, Jeffrey W., Hallett, Timothy B., & Garnett, Geoffrey P. (2011). Concurrent sexual partnerships and primary HIV infection: a critical interaction. Aids and behavior, 15(4), 687–692.
  • [\citenameFerguson & Garnett, 2000] Ferguson, Neil M., & Garnett, Geoffrey P. (2000). More realistic models of sexually transmitted disease transmission dynamics: sexual partnership networks, pair models, and moment closure. Sexually transmitted diseases, 27(10), 600–609.
  • [\citenameGoodreau et al., 2012] Goodreau, Steven M., Cassels, Susan, Kasprzyk, Danuta, Montaño, Daniel E., Greek, April, & Morris, Martina. (2012). Concurrent partnerships, acute infection and HIV epidemic dynamics among young adults in Zimbabwe. Aids and behavior, 16(2), 312–322.
  • [\citenameHallett et al., 2008] Hallett, Timothy B., Singh, Kanwarjit, Smith, Jennifer A., White, Richard G., Abu-Raddad, Laith J., & Garnett, Geoff P. (2008). Understanding the impact of male circumcision interventions on the spread of HIV in southern Africa. Plos one, 3(5), e2212.
  • [\citenameHamilton & Morris, 2010] Hamilton, Deven T., & Morris, Martina. (2010). Consistency of self-reported sexual behavior in surveys. Archives of sexual behavior, 39(4), 842–860.
  • [\citenameHamilton et al., 2008] Hamilton, Deven T, Handcock, Mark S, & Morris, Martina. (2008). Degree distributions in sexual networks: a framework for evaluating evidence. Sexually transmitted diseases, 35(1), 30.
  • [\citenameHolme & Saramäki, 2012] Holme, Petter, & Saramäki, Jari. (2012). Temporal networks. Physics reports, 519(3), 97–125.
  • [\citenameKarlin, 2014] Karlin, Samuel. (2014). A first course in stochastic processes. Academic press.
  • [\citenameKretzschmar & Dietz, 1998] Kretzschmar, Mirjam, & Dietz, Klaus. (1998). The effect of pair formation and variable infectivity on the spread of an infection without recovery. Mathematical biosciences, 148(1), 83–113.
  • [\citenameLeung & Kretzschmar, 2015] Leung, Ka Yin, & Kretzschmar, Mirjam. (2015). Concurrency can drive an HIV epidemic by moving R0 across the epidemic threshold:. Aids, 29(9), 1097–1103.
  • [\citenameLeung et al., 2014] Leung, Ka Yin, Kretzschmar, Mirjam, & Diekmann, Odo. (2014). SI infection on a dynamic partnership network: characterization of R_0. Journal of mathematical biology, 71(1), 1–56.
  • [\citenameLeung et al., 2012] Leung, K.Y., Kretzschmar, M.E.E., & Diekmann, O. (2012). Dynamic concurrent partnership networks incorporating demography. Theoretical population biology, 82(3), 229–239.
  • [\citenameMiller et al., 2011] Miller, Joel C., Slim, Anja C., & Volz, Erik M. (2011). Edge-based compartmental modelling for infectious disease spread. Journal of the royal society interface, rsif20110403.
  • [\citenameMoody, 2002] Moody, James. (2002). The importance of relationship timing for diffusion. Social forces, 81(1), 25–56.
  • [\citenameMorris, 1993] Morris, Martina. (1993). Telling tails explain the discrepancy in sexual partner reports. Nature.
  • [\citenameMorris & Kretzschmar, 1997] Morris, Martina, & Kretzschmar, Mirjam. (1997). Concurrent partnerships and the spread of HIV. Aids, 11(5), 641–648.
  • [\citenameMorris & Kretzschmar, 2000] Morris, Martina, & Kretzschmar, Mirjam. (2000). A microsimulation study of the effect of concurrent partnerships on the spread of HIV in Uganda. Mathematical population studies, 8(2), 109–133.
  • [\citenameMorris et al., 2009] Morris, Martina, Kurth, Ann E., Hamilton, Deven T., Moody, James, & Wakefield, Steve. (2009). Concurrent partnerships and HIV prevalence disparities by race: linking science and public health practice. American journal of public health, 99(6), 1023–1031.
  • [\citenameMorris et al., 2010] Morris, Martina, Epstein, Helen, & Wawer, Maria. (2010). Timing is everything: international variations in historical sexual partnership concurrency and HIV prevalence. Plos one, 5(11), e14092.
  • [\citenameNewman, 2007] Newman, M. E. J. (2007). Component sizes in networks with arbitrary degree distributions. Physical review e, 76(4), 045101.
  • [\citenameNicosia et al., 2012] Nicosia, Vincenzo, Tang, John, Musolesi, Mirco, Russo, Giovanni, Mascolo, Cecilia, & Latora, Vito. (2012). Components in time-varying graphs. Chaos: An interdisciplinary journal of nonlinear science, 22(2), 023101.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description