[
Abstract
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 closedform 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 crosssectional 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 crosssectional degree distribution (the number of partners at one instant) influences how connectivity emerges over time. Longterm 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 microlevel processes that govern node and link dynamics: the mean crosssectional degree , the edge dissolution rate , and the node exit rate . Conditional on these parameters, we vary the crosssectional 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 powerlaw distribution, are not represented in our analysis. However, longtailed distributions are usually found for cumulative degree (e.g., number of partners in the last year) [\citenameHamilton et al., 2008], rather than crosssectional degree. We will discuss the implications of longtailed 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 nonzero 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 pairformation 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 closedform 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 4tuple 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 crosssectional 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ősRé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 crosssectional network taken at equilibrium. Equivalently, can be the mean instantaneous degree of a node at equilibrium.
Given the set of node and edge activities, we can define a forward reachable path as a timeordered sequence of active edges between an initial and a destination node, within a time interval. An equivalent definition using crosssectional 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 crosssectional degree , crosssectional 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 nonzero 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 meanfield approximation (section 5). For the equilibrium persistence, we use the meanfield 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) 
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 crosssectional (static) network. For a Bernoulli degree distribution, clearly . For Poisson \shortcitenewman_component_2007, when
(1)  
(2)  
(3) 
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 nonhomogeneous 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 birthdeath 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 birthdeath 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 steadystate, , balancing the birth and death rates of the birthdeath process. For the Poisson case where all nodes are eligible for formation, this balance is . For the Bernoulli case, only the degree0 nodes are eligible to form an edge, and only degree1 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,
(4) 
in both the Poisson and Bernoulli cases and
(5) 
in the Poisson case and asymptotically in the Bernoulli case,
(6) 
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:
(7) 
due to standard properties of these distributions. In either the Poisson or Bernoulli case we can use iterated expectations to show that . Thus,
(8) 
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 [\citenameBenderdeMoll & 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 nonzero 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 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 takeoff, 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 vectorvalued 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 ,
(9) 
where we define
(10)  
(11) 
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
(12) 
where matrix
(13) 
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,
(14)  
(15) 
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:
(16)  
(17) 
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 degree0 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 degree0 node in the FRS. , , and at rate .
, . 
Death of a degree1 node in the FRS. , , and at rate .
, .
This leads to the following set of differential equations:
(18)  
(19)  
(20) 
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 crosssectional 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.
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 crosssectional 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.
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 .
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 meanfield models that modify the differential equations from Section 4.
For Poisson degree networks, the differential equation becomes
(21) 
For the Bernoulli case, we get
(22)  
(23) 
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 meanfield solutions for a 500 node network, up until the equilibrium.
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 meanfield equilibrium on the parameters.
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 meanfield model is deterministic and does not account for extinctions of the FRS, which happens stochastically. The meanfield solutions match up with the average persistent FRS size (for the set of initial nodes that do not have their FRS go extinct).
6.1 Persistent FRS size
The logistic meanfield equations from Section 5 allow us to calculate the equilibrium FRS fraction for the Poisson case, by setting the growth rate to zero:
(24) 
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:
(25) 
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 crosssectional 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
(26) 
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 meanfield models that . This is also illustrated in Figure 7 which shows the steadystate 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 degree0, which causes the network to have a nonhomogeneous degree distribution. The degree0 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.
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 nonFRS 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 nonnegative root of the equation [\citenameKarlin, 2014]
(27) 
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
(28) 
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 
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 nontrivial 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 nontrivial 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” (susceptibleinfected) 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; 3540 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 rightskewed exponential tail [\citenameLeung et al., 2012]. Mean crosssectional 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, powerlawlike 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 powerlaw network, we can get closer to a longtailed 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: groupspecific meandegree, 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 closedform results.
Additionally, a memoryless stochastic process (“Markov assumption”) does not account for different types of relationships. In real networks, partnership durations are often bimodal (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 nonFRS 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 noninfected. 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 timevarying 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], pairformation model [\citenameKretzschmar & Dietz, 1998], and edgebased 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 vectorvalued process,
(28)  
(29) 
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,
(30)  
(31) 
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 ,
(32) 
a.1 Poisson case
Applying this framework to the transitions in the Poisson case (see section 4.2), we derive the following equations:
(33)  
(34) 
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 closedform solutions
(35)  
(36) 
If we add in the equations relating to ,
(37)  
(38)  
(39)  
(40) 
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:
(41)  
(42)  
(43)  
(44)  
(45)  
(46)  
(47) 
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):
(48)  
(49)  
(50)  
(51)  
(52) 
Appendix B Adjustments to the Mean field Approximation
In our simulations, nodes enter the network with degree0, which causes the network to have a nonhomogeneous 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 degree0 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 nonFRS), the percent of singles is higher for the nonFRS 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  

“Single”  
“Paired”  
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  
NonFRS paired death 
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 meanfield 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 degree0 births is different for nodes in the FRS versus those not in the FRS. Since nodes enter the network as degree0 and nonFRS, the mean degree for nonFRS 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 meanfield equation (21) can be adjusted in the same way. We can obtain the mean component size for nonFRS 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.
References
 [\citenameBaggaley et al., 2006] Baggaley, Rebecca F, Garnett, Geoff P, & Ferguson, Neil M. (2006). Modelling the Impact of Antiretroviral Use in ResourcePoor 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.
 [\citenameBenderdeMoll & Morris, 2015] BenderdeMoll, 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 earlystage HIV infection does not predict the longterm 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., AbuRaddad, 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 selfreported 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). Edgebased 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 timevarying graphs. Chaos: An interdisciplinary journal of nonlinear science, 22(2), 023101.