The large graph limit of a stochastic epidemic model on a dynamic multilayer network ^{†}^{†}thanks: This research has been supported in part by the Mathematical Biosciences Institute and the National Science Foundation under grants RAPID DMS1513489 and DMS0931642.
Abstract
We consider an SIRtype (Susceptible Infected Recovered) stochastic epidemic process with multiple modes of transmission on a contact network. The network is given by a random graph following a multilayer configuration model where edges in different layers correspond to potentially infectious contacts of different types. We assume that the graph structure evolves in response to the epidemic via activation or deactivation of edges. We derive a large graph limit theorem that gives a system of ordinary differential equations (ODEs) describing the evolution of quantities of interest, such as the proportions of infected and susceptible vertices, as the number of nodes tends to infinity. Analysis of the limiting system elucidates how the coupling of edge activation and deactivation to infection status affects disease dynamics, as illustrated by a twolayer network example with edge types corresponding to community and healthcare contacts. Our theorem extends some earlier results deriving the deterministic limit of stochastic SIR processes on static, singlelayer configuration model graphs. We also describe precisely the conditions for equivalence between our limiting ODEs and the systems obtained via pair approximation, which are widely used in the epidemiological and ecological literature to approximate disease dynamics on networks. Potential applications include modeling Ebola dynamics in West Africa, which was the motivation for this study.
Keywords: Stochastic SIR process Configuration model Multilayer network Law of large numbers Ebola epidemic Multiple modes of transmission
MSC 92D30 60G55 60F
1 Introduction
A fundamental issue in disease dynamics is that contact patterns change in response to infection. This is particularly salient in the study of disease dynamics on contact networks: infected individuals curtail contacts with their regular community due to illness (e.g. being too sick to go to school or work) but increase their contacts with other segments of the population, such as healthcare workers or caretakers in the home. The recent Ebola outbreak in West Africa provides a stark example. The array and severity of symptoms, including high fever, diarrhea, vomiting, and hemorrhaging, make symptomatic individuals too ill to engage their regular community contacts and, instead, cause individuals to seek care in the home, hospital, or other facility. This coupling of evolution of network structure to disease status is basic, but a theoretical understanding of how this affects disease dynamics is lacking.
Disease dynamics on networks has been an extremely active area of research in the past 20 years, typically within the SIRtype (Susceptible Infected Recovered) modeling framework [14, 21, 32, 49, 56, 58, 74]. This has been stimulated in part by the explosion of data on networks of various sorts [5, 12, 13, 17, 24, 55, 65, 66, 76] and the recognition that network structure can have a dramatic impact on disease dynamics. Theoretical findings include applications of percolation theory to static networks [56, 30]. Less theory has been developed for networks that change over time, with much work in this area focusing on concurrent partnerships forming and breaking independent of disease status [1, 2, 40, 41, 72, 8, 23]. The study of adaptive networks, where the contact structure changes in response to disease progression, is an emerging area, as reviewed by Funk et al. [28]. One popular approach is to assume that susceptible individuals break connections to avoid infection [31, 25, 64, 79]. Related works examine behavioral changes due to awareness of infection [28, 29]. As these studies indicate, evolving network structure may lead to rich dynamics that are of practical importance for disease forecasting and evaluating public health interventions [63, 31, 64].
A challenge for understanding disease dynamics on networks is their high dimensionality – for example, a modestsized network or graph (here, and elsewhere in the paper, we use these terms interchangeably) may have tens of thousands of nodes and over a hundred thousand edges. Various approaches have been developed for deriving simpler models to approximate the full network dynamics including grouping vertices by degree [9, 10, 57] or “effective degree” [43, 45, 6] and considering stationary degree distributions for dynamic graphs [2, 40, 41]. Two approaches particularly relevant to this work are the pairwise approach (e.g. early work includes [59, 35, 36]) and the edgebased approach of Volz and Miller that is applicable to graphs with a specified degree distribution [71, 53, 54]. Both approaches naturally lead to consideration of the disease dynamics in the large graph limit, i.e. when the number of nodes tends to infinity. Whereas Volz and Miller derived their results heuristically, recent mathematical work has rigorously shown that a deterministic edgebased system of equations is the large graph limit of an SIR continuoustime Markov process on a static random graph [20, 33].
Multilayer networks, which allow for more complex disease dynamics, have also received much attention recently, as reviewed by Kivelä et al. [38]. In particular, multilayer networks where the interconnected layers can represent different populations have been considered [73, 60, 80, 15, 78]. The effect of degree correlation on twolayer networks has also been studied [62] with each layer an Erdős–Rényi or BarabásiAlbert random graph. Other models involving twolayer graphs were also considered where one layer corresponded to informationspreading and the second to disease transmission [34, 27, 26, 29] or where two competing pathogens spread on the two layers [75, 61]. Multilayer networks have also recently been employed to model temporal networks as sequences of static networks [67, 68]. The particular class of multilayer networks studied here are those in which the set of nodes is identical in each layer (i.e. nodealigned [38]) with distinct edge types corresponding to each layer. These are often referred to as multiplex (or multirelational) networks [19, 16, 77, 11] and can be represented as a graph with edges colored according to type.
In this work, we consider the problem of modeling an epidemic with different modes of disease transmission on a dynamic contact network. Specifically, we formulate an SIR continuoustime stochastic process on a multilayer graph, with specified degree distribution, where nodes represent individuals and edges represent potentially infectious contacts. Each layer contains the same set of nodes but corresponds to a different transmission mechanism (i.e. a multiplex network). In addition, we allow edges to be active or dormant with transmission occurring only along an active edge. The dynamic network structure allows us to incorporate behavioral changes due to infection. A simple example is a twolayer network with one layer corresponding to community contacts and the other to healthcare contacts where we assume that infected individuals deactivate their community edges, for example due to decreased mobility or isolation, while their healthcare edges are being activated as they seek care (Figure 1).
The main result of this work, Theorem 1 in Section 3, derives the large graph limit for the stochastic SIR process on the dynamic multilayer network. According to the theorem, the scaled counts of different edge and node types converge uniformly in probability to the solution of a deterministic system of equations. Thus, we obtain a relatively simple limiting model in the setting where network connectivity changes with the evolution of the disease process. In particular, it follows that for a certain class of random graphs the large graph limit coincides with the model obtained using either the pairwise or edgebased approach. As we demonstrate with the twolayer network example, the limiting systems are amenable for mathematical analysis, allowing us to gain biological insight into how changing network structure influences disease dynamics. Moreover, Theorem 1 extends previous results on edgebased models [20, 33].
This paper is organized as follows. Section 2 introduces the stochastic model that is considered along with the necessary notation. Section 3 presents our main result, a law of large numbers for the stochastic process on the dynamic multilayer network, and considers several important special cases, which relate our result to edgebased and pairwise models. The twolayer communityhealthcare network model and its analysis is given in Section 4. The proof of our main theorem is given in Section 5. We conclude with a discussion in Section 6. The Appendices provide further mathematical details.
2 Stochastic model
We start by introducing some notation and relevant definitions for dynamic multilayer networks. This includes defining the class of random graphs to be considered and extending the notions of degree and excess degree to the multilayer setting. Then, we introduce the stochastic process considered in this work, which is the appropriately modified version of the SIR process.
2.1 Layered configuration model
Let denote the number of layers and, for any vectors in , denote . The probability generating function (pgf) of the multivariate degree distribution is given by
(1) 
where is the probability of a node being of (multi)degree , i.e. having neighbors in layer .
Given a realization of the degree distribution on nodes, we construct a multilayer graph as follows. Each node is assigned a collection of halfedges in each layer corresponding to its degree, and then halfedges within each layer are paired uniformly at random. Thus, restricted to the th layer, the resulting graph is a realization of a configuration model [55] with the degree distribution given by the th marginal of . We refer to the collection of such realized graphs as the layered configuration model (LCM) and denote it by .
Excess degree distributions.
In a singlelayer graph, the excess degree of a node is calculated by following an edge to from a neighbor and counting the number of other neighbors (not including ) of [55]. It will be convenient to extend the notion of an excess degree distribution to the multilayer setting. Let denote the probability that a randomly selected neighbor (i.e. neighbor in layer ) of a node has degree (i.e. degree in layer ) equal to . Then, by LCM construction, is given as
where is the average degree, denotes the partial derivative with respect to , and is the vector of ones in . Correspondingly, let denote the pgf of the excess degree distribution of a node randomly selected as a neighbor. Then,
where is the vector of ones with the th coordinate replaced by . The average excess degree of an neighbor is then given by
Finally, we define the normalized average excess degree of an neighbor as
(2) 
Note that, for the univariate (i.e. single layer) case when , , which is the wellknown distribution of the degree of a neighbor (also referred to as the sizebiased distribution [70] and corresponding to the excess degree distribution [55]), and is the ratio of the average excess degree to average degree.
2.2 process
Assume that we have a realization of an LCM specifying the contact network for a population of size . The disease modeling framework adopted is the standard SIR compartmental model where individuals are classified based on their infection status [37]. , and correspond, respectively, to susceptible, infected, and recovered (or removed) individuals. We assume that edges within layers represent potentially infectious contacts of a certain type and we allow the network to be dynamic in response to infection. That is, we assume that nodes will either activate or deactivate their edges, depending on their infection status. An infectious node drops (resp. activates) edges in layer at rate (resp. ). We assume that a layer cannot be both activating and deactivating, i.e. at most one of and are nonzero, and we also assume that all deactivating layer edges are initially activated and all activating layer edges are initially deactivated. Let denote the deactivating layers (with ) and let denote the activating layers (with ). Then, event types may occur: infection (I) along an edge of any of the types, drop () of a deactivating edge or activation () of an activating edge, and recovery (). The timings of all events are assumed to follow independent exponential clocks with the following rates:
rate of infection along  edges (),  
rate of deactivation (drop) of  edges,  
rate of activation of  edges,  
rate of recovery (). 
For a susceptible node , let and denote, respectively, the number of infectious and susceptible active neighbors of . Similarly, let and denote the number of deactivated neighbors of . Also, for an infected node , let and denote the number of susceptible active and deactivated, respectively, neighbors of . We consider aggregate variables that are the total number of nodes or pairs of neighboring nodes (i.e. dyads) with a given disease status. For example, the total number of edges between susceptible and infectious individuals is denoted and is given by . We denote the aggregate dyad counts as vectors in , e.g. and likewise for , , and . Note that and count the edges twice. We let denote the state of the aggregate stochastic process at time where and denote the number of susceptible and infectious nodes, respectively. Note that the number of recovered individuals is given by and so, for the sake of simplicity, we ignore the equation for throughout. The transitions for the aggregate process are listed in Table 1. We refer to such a process as in order to emphasize the activation and deactivation events. The analysis of this process is complicated, partially due to the aggregation of the nodes that destroys the Markov property (see, e.g. [2]).
Event  Rate  Transition  Arrangement 
Infection of by along edge,  
for  
for  
for  
for  
Deactivation of edge,  
Activation of edge,  
Recovery of infected  
for  
for 
Note that, as above, the notation of a dyad subscript (e.g. ) is understood throughout to denote a (row) vector in . Also, with a slight abuse of notation we take multiplication, division, integration and ordering of vectors to be coordinatewise. The state variables depend on but we do not explicitly acknowledge this in our notation.
Consider the process on the LCM with transitions as outlined in Table 1. The DoobMeyer decomposition theorem [48] guarantees the existence of a zeromean martingale such that
(3) 
where the integrable function is given by
(4)  
We now define two more variables that will help us describe the evolution of the process in the large graph limit. Let where is the number of edges belonging to susceptible nodes at time . We partition the collection of susceptible nodes by their degree so that which corresponds also to Note
(5) 
We also define by
(6) 
where . We may interpret as the probability of no infection along a edge by time in a susceptible node of degree one, given that the node was susceptible at . That is, is the probability that a susceptible node of (multi)degree has not been infected through any layer by time , given that the node was susceptible at . Note also that we may equivalently write
(7) 
where and
(8) 
As shown in Theorem 1 below, plays a key role in describing the evolution of in the large graph limit. The use of such a variable was pioneered by Volz [71] and Miller [50] in their edgebased approach. In fact, as shown in Section 3.2, in the singlelayer, static network case, the large graph limit of corresponds to the variable in the standard SIR edgebased model [53].
3 Large graph limit theorems
The stochastic process defined in Section 2.2 is complex and difficult to analyze. In this section we present a limit theorem (Theorem 1) that shows, under mild technical assumptions, that the stochastic process converges to a relatively simple system of ODEs as the number of nodes tends to infinity. The limiting ODEs retain key features of the epidemic process while being amenable to analysis. In the case of a finite but large population, analysis of this deterministic system provides a good approximation to disease dynamics. In the general case, the evolution of the quantities of interest, , will involve a function of the variable defined in the preceding section. In Section 3.2, we state corollaries that relate our result to edgebased models in the special case of static graphs [54]. Finally, in Section 3.3 we show that, for a particular class of degree distributions, the evolution of decouples from , revealing a perhaps surprising connection between our limiting system and lowdimensional models obtained via pair approximation.
3.1 General case
All limits considered below, unless otherwise noted, are for . We use to denote convergence in probability in the space of rightcontinuous with finite left limits (càdlàg) stochastic processes on random configurations drawn according to an LCM . We say that a sequence of random variables with high probability (w.h.p.) if for any . Let . We make the following assumptions:

For ,

The fractions of initially susceptible, infected, and recovered nodes converge, respectively, to some , , , i.e.
Furthermore, , , and the initially infected and recovered nodes are chosen randomly.

The assumption 1 implies that, for large graphs, the infection does not overtake the graph. That is, some proportion of the individuals remain susceptible on and, hence, is welldefined in (6). Furthermore, it implies that the average degree of a randomly chosen node, i.e. , is positive since .
Assumption 2 implies that the initial conditions for the dyads, scaled by , also converge in probability, i.e.
(9) 
where, for the deactivating layers ,
and, for the activating layers ,
For example, we obtain the initial condition for as follows. By assumption, all deactivating layer edges are initially activated and all activating layer edges are initially deactivated. Then, according to 2, the limiting probability of selecting a random node that is susceptible is . The average number of edges a node has is , and the limiting probability that a given edge connects to an infected node is . Therefore, . The other dyad initial conditions are obtained similarly. We denote in what follows.
The assumption 3 implies that for (i.e. the second moments of the marginal degree distributions are finite) and also that the multigraph considered when matching halfedges uniformly at random is a simple graph with positive probability [33].
Before stating the main result, we define a quantity that plays a key role in describing how the network structure affects the large graph limit. For , let be defined by
(10) 
As we discuss in Section 3.3.3, can be interpreted as the ratio of the average excess degree of a susceptible node chosen randomly as an neighbor of an infectious node to the average degree of a susceptible node.
We also define the function that, by Theorem 1, describes the evolution of in the large graph limit. Let and define where and are given by
(11)  
Theorem 1 (Strong law of large numbers).
3.2 Edgebased limiting systems
We consider two special cases of the large graph limit theorem for multilayer networks. First, we consider a static network, i.e. the case where for . Corollary 1 shows that in this case our system (12) is equivalent to an edgebased model with multiple modes of transmission. The model is one proposed by Miller and Volz [54] but with a modification to allow for a large number of initially infected nodes (following [51] where Miller modifies the standard SIR edgebased model for such a scenario). In the case that the initially infected nodes are randomly chosen (which we assume in 2), the model is given by
(13) 
Details of the equivalency, i.e. proof of Corollary 1, are given in Appendix A.
Corollary 1.
We further consider the special case of a static, single layer graph (i.e. and ). In the case , (13) reduces to the wellknown edgebased SIR model of Volz and Miller et al. [53, 71], which has been proven to be the large graph limit of the SIR stochastic process on a static configuration model graph [33, 20]. By taking in Corollary 1 we have provided an alternative proof of this fact.
3.3 Pairwise limiting systems
In this section we consider a certain class of LCMs where as defined by equation (10) is constant. This affords a substantial simplification to the limiting system (11) and, in fact, the system of differential equations defining the large graph limit will coincide with the model derived via the pairwise approach.
3.3.1 Poissontype distributions
We define a multivariate Poissontype (PT) distribution to be a distribution with a pgf that satisfies
(14) 
where we recall the definition of the normalized average excess degree, , in equation (2). At first glance, (14) may seem opaque; however, in fact it is satisfied by a broad class of distributions. For example, in the single layer case this condition is equivalent to
which is satisfied by the univariate Poisson (), binomial (), and negative binomial () distributions. Note that a regular graph (where all nodes have degree ) is a special case of the binomial distribution. Additionally, the geometric distribution is a special case of the negative binomial distribution. Bansal et al. have shown that the geometric distribution (i.e. the discrete analog of the exponential distribution) gives the best fit for several empirical contact networks [7]. In the multilayer case with , if the layers are independent, i.e. the pgf can be written as , then the PT condition (14) reduces to each layer being a (univariate) Poissontype distribution.
3.3.2 Pairwise model
If the degree distribution is PT as defined by condition (14), the limiting system (11) defining has a particularly simple form. Indeed, substituting the constant for decouples from . We consider the resulting model in this section and introduce some new notation to do so. Let and denote, respectively, the number of activated and deactivated edges of type between a node of status and a node of status . Let denote the number of triples with an active edge between nodes of status and and an active edge between node the node of status and one of status . Similarly, will denote such triples where the edge is dormant.
Under the correlation equations approach of Rand [59], triples are needed to describe the evolution of pairs, quadruples (e.g. ) are needed to describe the evolution of triples, and so forth. A pair approximation for triples is used in order to close the system at the level of pairs [59]. For consideration of triples in the multilayer setting, we must take into account the edge types and the appropriate excess degrees as defined in Section 2.1. Let denote the probability that a node has disease status given an edge (or triple) arrangement for . The pair approximation of is then calculated as follows:
(15) 
with as defined in equation (2). We can alternatively calculate to show that .
Applying the correlation equations approach using the pair approximation (15) to the dynamics described in Section 2.2 results exactly in the same equations as the limiting system (12) in the case of a PT distribution:
(16)  
Here, we have derived the system (16) using absolute pair and triple counts. Notice that, if we scale all variables by the graph size , the nondimensional quantities satisfy the same system of equations (this holds for any and, hence by continuity of the solution, also in the limit). From here on we will consider only the scaled variables, which will be convenient when we state the law of large numbers in Corollary 3. Accordingly, we set the initial conditions to be
(17) 
so that they agree with the initial conditions in Theorem 1 as .
Observe also that, in fact, we can reduce the dimension of the system (16) since we only need to keep track of the deactivated edges for the activating layers, i.e. . Also, for an activating layer since its initial condition is zero and, hence, we only must track for . We refer to system (16) with its initial condition (17) as the pairwise model.
The discussion above is summarized in the following corollary.
Corollary 3.
Furthermore, we can consider the implications of Corollary 3 in the static, single layer case. If and , then the pairwise model (16) reduces to the wellknown correlation equations model of Rand [59]:
(18)  
Corollary 4.
Note that together Corollaries 1 and 3 imply that, in the case of a multivariate PT degree distribution on a static graph, the pairwise model (16) is equivalent to the edgebased model with multiple modes of transmission (13). Likewise, Corollaries 2 and 4 indicate that the pairwise SIR model (18) is equivalent to the edgebased SIR model, (13) with , when the distribution is PT. We note that the edgebased SIR model has previously been shown to be equivalent [32, 52] to a higher dimensional pairwise model of Eames and Keeling [22]. The latter model stratifies the susceptible nodes by degree and, hence, has dimension where represents the number of distinct degrees [52]. The model of dimension was derived as an approximation to an earlier wellknown model of Eames and Keeling, which is of dimension [22]. We observe that, in fact, (18) can be reduced to two differential equations. Separation of variables on (see, e.g., [3] for ) gives . Subsequent inspection of