The large graph limit of a stochastic epidemic model on a dynamic multilayer network

ABSTRACT

We consider a Markovian SIR-type (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 of infectious nodes. 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 two-layer network example with edge types corresponding to community and healthcare contacts. Our theorem extends some earlier results describing the deterministic limit of stochastic SIR processes on static, single-layer 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. The flexible modeling framework and asymptotic results have potential application to many disease settings including Ebola dynamics in West Africa, which was the original motivation for this study.

KEYWORDS:

• Stochastic SIR process
• configuration model
• multilayer network
• law of large numbers
• Ebola epidemic
• multiple modes of transmission

AMS SUBJECT CLASSIFICATIONS:

• 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 currently lacking.

Disease dynamics on networks has been an extremely active area of research in the past 20 years, typically within the SIR-type (Susceptible → Infected → Recovered) modeling framework [17,25,42,65,71,75,100]. This has been stimulated in part by the explosion of data on networks of various sorts [6,14,16,21,30,72,89,91,105] and the recognition that network structure can have a dramatic impact on disease dynamics. Theoretical findings include applications of percolation theory to static networks [39,71]. 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,9,29,54,55,98]. 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. [33]. One popular approach is to assume that susceptible individuals break connections to avoid infection [31,40,88,109]. Related works examine behavioural changes due to awareness of infection [33,38]. 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 [40,85,88].

A challenge for understanding disease dynamics on networks is their high dimensionality – for example, a modest-sized 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 [11,12,74] or ‘effective degree’ [7,57,59] and considering stationary degree distributions for dynamic graphs [2,54,55]. Two approaches particularly relevant to this work are the pairwise approach (e.g. early work includes [45,46,78]) and the edge-based approach of Volz and Miller that is applicable to graphs with a specified degree distribution [69,70,97]. 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 edge-based system of equations is the large graph limit of an SIR continuous-time Markov process on a static random graph [24,43].

Multilayer networks, which allow for more complex disease dynamics, have also received much attention recently, as reviewed by Kivelä et al. [48]. In particular, multilayer networks where the interconnected layers can represent different populations have been considered [18,81,99,108,110]. The effect of degree correlation on two-layer networks has also been studied [84] with each layer being an Erd$ő$os–Rényi or Barabási-Albert random graph. Other models involving two-layer graphs were also considered where one layer corresponded to information-spreading and the second to disease transmission [34,35,38,44] or where two competing pathogens spread on the two layers [82,102]. Multilayer networks have also recently been employed to model temporal networks as sequences of static networks [92,93]. The particular class of multilayer networks studied here are those in which the set of nodes is identical in each layer (i.e. node-aligned [48]) with distinct edge types corresponding to each layer. These are often referred to as multiplex (or multi-relational) networks [13,19,23,107] 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 continuous-time 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 network structure is dynamic in that edges can activate or deactivate over the course of infection. This approach allows us to incorporate behavioural changes due to infection while keeping the total edges (active plus dormant) given by the degree distribution fixed. A simple example is a two-layer 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).

Figure 1. (Left) Neighborhood of a susceptible vertex (labeled 1) with an infected (red) neighbour. Community (top/blue) and healthcare (bottom/green) contacts are shown as active (solid) or deactivated (dashed). (Right) After infection of vertex 1, two of its healthcare contacts are activated and one community contact is dropped.

The main result of this work, Theorem 3.1 in Section 3, describes 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 [45] or the edge-based [98] approximation approach. As we demonstrate with the two-layer 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 3.1 extends previous results on edge-based models [24,43].

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 edge-based and pairwise models. The two-layer community-healthcare network model and its analysis is given in Section 4. We conclude with a discussion in Section 5. The proof of our main theorem is given in Appendix 2 which also provides further mathematical details along with a summary of notation for the main body of the paper.

2. Stochastic model

Recent advances in computational methods and the ever-increasing power of modern computers have made it possible to consider stochastic versions of the classical ODE epidemic models. Such models not only provide the overall trend of an epidemic across a population but also inform about the stochastic fluctuations around the mean and, hence, about the intrinsic noise in the system (see, for instance [106] and references therein). Some visual examples are provided in the next section. The stochastic models are typically formulated as continuous time Markov processes with discrete state space (see, e.g. [1,43,97]), which is also the framework we adopt in this paper.

We start by introducing some notation and relevant definitions for dynamic multilayer networks. The class of random graphs considered here will be an extension of the configuration model to the multilayer setting. In the following section we will precisely define the layered configuration model, and extend 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 r denote the number of layers and, for any vectors $\mathit{x}=(x_1,\dots ,x_r),\mathit{k}=(k_1,\dots ,k_r)$ in ${\mathbb{R}}^r$, denote ${\mathit{x}}^{\mathit{k}}=\prod _{i=1}^r{x}_i^{k_i}$. The probability generating function (pgf) of the multivariate degree distribution is given by (1) $\psi \left(\mathit{x}\right)=\sum _{\mathit{k}}{p}_{\mathit{k}}{\mathit{x}}^{\mathit{k}}$(1) where $p_{\mathit{k}}=P(k_1,\dots ,k_r)$ is the probability of a node being of (multi-)degree $\mathit{k}$, i.e. having $k_i$ neighbours in layer i.

Given a realization of the degree distribution on n nodes, we construct a multilayer graph as follows. Each node is assigned a collection of half-edges in each layer corresponding to its degree, and then half-edges within each layer are paired uniformly at random. We assume that the pairing is done independently in each layer. Thus, restricted to the jth layer, the resulting graph is a realization of a configuration model (see, Chapter 13 in [72]) with the degree distribution given by the jth marginal of ψ (also see, for instance, Section 2.2.4 in [69]). We refer to the collection of such realized graphs as the layered configuration model (LCM) and denote it by ${\mathcal{G}}_r(\psi ,n)$.

Excess degree distribution

In a single-layer graph, the excess degree of a node u is calculated by following an edge to u from a neighbour v and counting the number of other neighbours (not including v) of u (see, Chapter 13 in [72]). It will be convenient to extend the notion of an excess degree distribution to the multilayer setting. Let $P_{j|i}\left(l\right)$ denote the probability that a randomly selected i-neighbour (i.e. neighbour in layer i) of a node u has j-degree (i.e. degree in layer j) equal to l. Then, by LCM construction, $P_{j|i}\left(l\right)$ is given as $P_{j|i}\left(l\right)=\sum _{\mathit{k}:k_j=l}{k}_i{p}_{\mathit{k}}/{\mu }_i$ where ${\mu }_i={\mathrm{\partial }}_i\psi \left(\mathbf{1}\right)=\sum _{\mathit{k}}{k}_i{p}_{\mathit{k}}$ is the average i-degree, ${\mathrm{\partial }}_i$ denotes the partial derivative with respect to $x_i$, and $\mathbf{1}$ is the vector of ones in ${\mathbb{R}}^r$. Correspondingly, let ${\psi }_{j|i}^{ex}$ denote the pgf of the excess j-degree distribution of a node randomly selected as a i-neighbour. Then, ${\psi }_{j|i}^{ex}\left(x_j\right)=\sum _l{P}_{j|i}\left(l\right)x_j^l=\sum _{\mathit{k}}\frac{k_i{p}_{\mathit{k}}}{{\mu }_i}{x}_j^{k_j}=\frac{1}{{\mu }_i}{\mathrm{\partial }}_i\psi (1,\dots ,1,x_j,1,\dots ,1)$ where $(1,\dots ,1,x_j,1,\dots ,1)$ is the vector of ones with the jth coordinate replaced by $x_j$. The average excess j-degree of an i-neighbour is then given by ${\mu }_{j|i}^{ex}={\mathrm{\partial }}_j{\psi }_{j|i}^{ex}\left(1\right)=\frac{1}{{\mu }_i}{\mathrm{\partial }}_{ij}^2\psi \left(\mathbf{1}\right).$ Finally, we define the normalized average excess j-degree of an i-neighbour as (2) ${\kappa }_{ji}=\frac{{\mu }_{j|i}^{ex}}{{\mu }_j}=\frac{{\mathrm{\partial }}_{ij}^2\psi \left(\mathbf{1}\right)}{{\mathrm{\partial }}_i\psi \left(\mathbf{1}\right){\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}.$(2) Note that, for the univariate (i.e. single layer) case when r=1, $P_{i|i}\left(k\right)=k{p}_k/\mu$, which is the well-known distribution of the degree of a neighbour (also referred to as the size-biased distribution [96] and corresponding to the excess degree distribution $q_k=(k+1)p_{k+1}/\mu$ [72]), and κ is the ratio of the average excess degree to average degree.

2.2. SIdaR process

Assume that we have a realization of an LCM ${\mathcal{G}}_r(\psi ,n)$ specifying the contact network for a population of size n. The disease modeling framework adopted is the standard Markovian SIR compartmental model where individuals are classified based on their infection status [47]. S, I and R 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 infected nodes will either activate or deactivate their edges, depending on edge type. An infectious node drops (resp. activates) edges in layer j at rate ${\delta }_j$ (resp. ${\eta }_j$). We assume that a layer cannot be both activating and deactivating, i.e. at most one of ${\delta }_j$ and ${\eta }_j$ are nonzero, and we also assume that all deactivating layer edges are initially activated and all activating layer edges are initially deactivated. Note that both active and deactivated edges are counted in a node's degree, which is therefore constant throughout the course of an epidemic (see Section 2.3.1 in [69] for a similar approach). Let $j=1,\dots ,m$ denote the deactivating layers (with ${\eta }_j=0$) and let $j=m+1,\dots ,r$ denote the activating layers (with ${\delta }_j=0$). Then, 2r+1 event types may occur: infection (I) along an edge of any of the r types, drop (d) of a deactivating edge or activation (a) of an activating edge, and recovery (R). The timings of all events are assumed to follow independent exponential clocks with the following rates: $\begin{array}{lll}{\beta }_j& \text{rate of infection along }j\text{-edges }\left(S\stackrel{j}{⟶}I\right),& j=1,\dots ,r\\ {\delta }_j& \text{rate of deactivation (drop) of }j\text{-edges},& j=1,\dots ,m\\ {\eta }_j& \text{rate of activation of }j\text{-edges},& j=m+1,\dots ,r\\ \gamma & \text{rate of recovery }(I⟶R).& \end{array}$

For a susceptible node u, let $X_j^{SI,u}$ and $X_j^{SS,u}$ denote, respectively, the number of infectious (i.e. infected, not yet recovered) and susceptible active j-neighbours of u. Similarly, let $X_j^{\widetilde{SI},u}$ and $X_j^{\widetilde{SS},u}$ denote the number of deactivated j-neighbours of u. Also, for an infected node u, let $X_j^{IS,u}$ and $X_j^{\widetilde{IS},u}$ denote the number of susceptible active and deactivated, respectively, j-neighbours of u. We consider aggregate variables that are the total number of nodes or pairs of neighbouring nodes (i.e. dyads) with a given disease status. For example, the total number of j-edges between susceptible and infectious individuals is denoted $X_j^{SI}$ and is given by $X_j^{SI}=\sum _{u\in S}{X}_j^{SI,u}$. We regard the aggregate dyad counts as vectors in ${\mathbb{R}}^r$, e.g. ${\mathit{X}}^{SI}=(X_1^{SI},\dots ,X_r^{SI})$ and likewise for ${\mathit{X}}^{\widetilde{SI}}$, ${\mathit{X}}^{SS}$, and ${\mathit{X}}^{\widetilde{SS}}$. Note that ${\mathit{X}}^{SS}$ and ${\mathit{X}}^{\widetilde{SS}}$ count the edges twice. We let $\mathit{X}\left(t\right)=(X^S,X^I,{\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SI}},{\mathit{X}}^{SS},{\mathit{X}}^{\widetilde{SS}})\left(t\right)$ denote the state of the aggregate stochastic process at time t>0 where $X^S$ and $X^I$ denote the number of susceptible and infectious nodes, respectively. Note that the number of recovered individuals is given by $X^R=n-X^S-X^I$ and so, for the sake of simplicity, we ignore the equation for $X^R$ throughout. In addition, we do not keep track of ${\mathit{X}}^{\mathit{I}\mathit{I}}$ or the dyads of recovered individuals since the evolutions of the main quantities of interest, S and I, are not affected by these variables. The transitions for the aggregate process are listed in Table 1. We refer to such a process as $SIdaR(r,m)$ 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]).

Table 1. Transitions for the $SIdaR(r,m)$ process according to the 2r+1 possible event types with corresponding rates. Network arrangements corresponding to the transitions are also given with $\stackrel{j}{-}$ and $\stackrel{j}{\sim }$ denoting, respectively, active and deactivated edges of type j between nodes (denoted u, v and w). Here, $N_l^{Y,u}$ denotes the set of l-neighbours of node u with disease status Y .

Event	Rate	Transition	Arrangement
Infection of u by v	${\beta }_j{X}_j^{SI}$	$(X^S,X^I)\to (X^S-1,X^I+1)$
along j-edge		$(X_l^{SS},X_l^{SI})\to (X_l^{SS}-X_l^{SS,u},X_l^{SI}+X_l^{SS,u})$	$v\stackrel{j}{-}u\stackrel{l}{-}w$  for $w\in N_l^{S,u}$
$j=1,\dots ,r$		$(X_l^{\widetilde{SS}},X_l^{\widetilde{SI}})\to (X_l^{\widetilde{SS}}-X_l^{\widetilde{SS},u},X_l^{\widetilde{SI}}+X_l^{\widetilde{SS},u})$	$v\stackrel{j}{-}u\stackrel{l}{\sim }w$ for $w\in N_l^{S,u}$
		$X_l^{SI}\to X_l^{SI}-X_l^{SI,u}$	$v\stackrel{j}{-}u\stackrel{l}{-}w$ for $w\in N_l^{I,u}$
		$X_l^{\widetilde{SI}}\to X_l^{\widetilde{SI}}-X_l^{\widetilde{SI},u}$	$v\stackrel{j}{-}u\stackrel{l}{\sim }w$ for $w\in N_l^{I,u}$
Deactivation of j-edge	${\delta }_j{X}_j^{SI}$	$X_j^{SI}\to X_j^{SI}-1$
$j=1,\dots ,m$
Activation of j-edge	${\eta }_j{X}_j^{\widetilde{SI}}$	$(X_j^{\widetilde{SI}},X_j^{SI})\to (X_j^{\widetilde{SI}}-1,X_j^{SI}+1)$
$j=m+1,\dots ,r$
Recovery of infected u	$\gamma X^I$	$X^I\to X^I-1$
		$X_j^{SI}\to X_j^{SI}-X_j^{IS,u}$	$u\stackrel{j}{-}w$  for $w\in N_l^{S,u}$
		$X_j^{\widetilde{SI}}\to X_j^{\widetilde{SI}}-X_j^{\widetilde{IS},u}$	$u\stackrel{j}{\sim }w$  for $w\in N_l^{S,u}$

Note that the dyad (e.g. ${\mathit{X}}^{SI}$) is understood throughout as a (row) vector in ${\mathbb{R}}^r$. Also, with a slight abuse of notation we take multiplication, division, integration and ordering of vectors to be coordinatewise. The state variables depend on n but we do not explicitly acknowledge this in our notation.

Consider the $SIdaR(r,m)$ process $\mathit{X}\left(t\right)$ on the LCM ${\mathcal{G}}_r(\psi ,n)$ with transitions as outlined in Table 1. The Doob-Meyer decomposition theorem [64] guarantees the existence of a zero-mean martingale $\mathit{M}\left(t\right)=(M^S,M^I,{\mathit{M}}^{SI},{\mathit{M}}^{\widetilde{SI}},{\mathit{M}}^{SS},{\mathit{M}}^{\widetilde{SS}})\left(t\right)$ such that (3) $\mathit{X}\left(t\right)=\mathit{X}\left(0\right)+{\int }_0^t{\mathit{F}}^{\mathit{X}}\left(\mathit{X}\right(s\left)\right)\mathrm{d}s+\mathit{M}\left(t\right)$(3) where the integrable function ${\mathit{F}}^{\mathit{X}}\left(\mathit{X}\right)=({\mathcal{F}}^S,{\mathcal{F}}^I,{\mathit{F}}^{\mathit{S}\mathit{I}},{\mathit{F}}^{\widetilde{\mathit{S}\mathit{I}}},{\mathit{F}}^{\mathit{S}\mathit{S}},{\mathit{F}}^{\widetilde{\mathit{S}\mathit{S}}})\left(\mathit{X}\right)$ is given by (4) $\begin{array}{rlr}{\mathcal{F}}^S\left({\mathit{X}}^{SI}\right)& =-\sum _{l=1}^r{\beta }_l{X}_l^{SI},& \\ {\mathcal{F}}^I(X^I,{\mathit{X}}^{SI})& =\sum _{l=1}^r\left({\beta }_l{X}_l^{SI}\right)-\gamma X^I,& \\ {\mathcal{F}}_j^{SI}({\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SI}},{\mathit{X}}^{SS})& =\sum _{i\in S}\left(\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}(X_j^{SS,i}-X_j^{SI,i})\right)-(\gamma +{\delta }_j)X_j^{SI}+{\eta }_j{X}_j^{\widetilde{SI}},\\ {\mathcal{F}}_j^{\widetilde{SI}}({\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SI}},{\mathit{X}}^{\widetilde{SS}})& =\sum _{i\in S}\left(\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}(X_j^{\widetilde{SS},i}-X_j^{\widetilde{SI},i})\right)-(\gamma +{\eta }_j)X_j^{\widetilde{SI}}+{\delta }_j{X}_j^{SI},\\ {\mathcal{F}}_j^{SS}({\mathit{X}}^{SI},{\mathit{X}}^{SS})& =-2\sum _{i\in S}\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}{X}_j^{SS,i},\\ {\mathcal{F}}_j^{\widetilde{SS}}({\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SS}})& =-2\sum _{i\in S}\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}{X}_j^{\widetilde{SS},i},\end{array}$(4) for $j=1,\dots ,r$.

We now define two more variables that will help us describe the evolution of the process in the large graph limit. Let ${\mathit{X}}^{S\text{•}}\left(t\right)=(X_1^{S\text{•}},\dots ,X_r^{S\text{•}})\left(t\right)$ where $X_j^{S\text{•}}\left(t\right)$ is the number of j-edges (active and deactivated) belonging to susceptible nodes at time t. We partition the collection of susceptible nodes S by their degree $\mathit{k}\ge 0$ so that $S=\cup S_{\mathit{k}}$, which corresponds also to $X^S=\sum _{\mathit{k}}{X}^{S_{\mathit{k}}}.$ Note (5) ${\mathit{X}}^{S\text{•}}\left(t\right)=\sum _{\mathit{k}}\mathit{k}{X}^{S_{\mathit{k}}}\left(t\right).$(5) We also define $\mathbf{\Theta }\left(t\right)=({\mathrm{\Theta }}_1,\dots ,{\mathrm{\Theta }}_r)\left(t\right)$ by (6) $\mathbf{\Theta }\left(t\right)=\exp \left(-\mathit{\beta }{\int }_0^t\frac{{\mathit{X}}^{SI}\left(s\right)}{{\mathit{X}}^{S\text{•}}\left(s\right)}\mathrm{d}s\right)$(6) where $\mathit{\beta }=({\beta }_1,\dots ,{\beta }_r)$. We may interpret ${\mathrm{\Theta }}_j\left(t\right)$ as the probability of no infection along a j-edge by time t in a susceptible node of j-degree one, given that the node was susceptible at t=0. That is, ${\mathbf{\Theta }}^{\mathbf{1}}=\prod _{j=1}^r{\mathrm{\Theta }}_j$ is the probability that a susceptible node of (multi-)degree $\mathbf{1}$ has not been infected through any layer by time t, given that the node was susceptible at t=0. Note also that we may equivalently write (7) $\mathbf{\Theta }\left(t\right)=\mathbf{\Theta }\left(0\right)+{\int }_0^t{\mathit{F}}^{\mathbf{\Theta }}\left({\mathit{X}}^{SI}\right(s),{\mathit{X}}^{S\text{•}}(s),\mathbf{\Theta }(s\left)\right)\mathrm{d}s$(7) where $\mathbf{\Theta }\left(0\right)=\mathbf{1}$ and (8) ${\mathit{F}}^{\mathbf{\Theta }}({\mathit{X}}^{SI},{\mathit{X}}^{S\text{•}},\mathbf{\Theta })=-\mathit{\beta }\mathbf{\Theta }{\mathit{X}}^{SI}/{\mathit{X}}^{S\text{•}}.$(8) As shown in Theorem 3.1 below, $\mathbf{\Theta }\left(t\right)$ plays a key role in describing the evolution of $\mathit{X}\left(t\right)$ in the large graph limit. The use of such a variable was pioneered by Volz [97] and Miller [66] in their edge-based approach. In fact, as shown in Section 3.2, in the single-layer, static network case, the large graph limit of $\mathbf{\Theta }\left(t\right)$ corresponds to the variable in the standard SIR edge-based model [70].

3. Large graph limit theorems

The stochastic process defined in Section 2.2 is complex and difficult to analyze directly. In this section we present a limit theorem (Theorem 3.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 sense that fluctuations around the mean due to intrinsic noise in the system shrink as the graph size grows. The study of such large volume limits for stochastic processes of the type discussed here was originated by Kurtz in [50] in the context of chemical reaction models. His work has subsequently inspired multiple large volume results on the stochastic SIR-type models (see, e.g. Chapter 5 in [4]). Among others, Andersson [3] derived limit theorems for a discrete-time random graph epidemic model under rather restrictive assumptions on the degree sequence of the random graph, such as finiteness of a $(4+ϵ)$th moment for some $ϵ>0$. Using a heuristic argument, Volz [97] presented scaling limits for an SIR model on random graphs in the form of ODEs. Decreusefond et al. [24] later proved Volz's results rigorously by summarizing the epidemic dynamics on a configuration model random graph using certain measure-valued process. Several similar law of large numbers-type scaling limits under varying sets of technical assumptions surfaced afterwards. For example, Bohman and Picollelli [15] and Barbour and Reinert [10] assumed uniformly bounded degrees. Janson, Luczak and Windridge [43] assumed the degree of a randomly chosen susceptible vertex to be uniformly integrable and the maximum degree of initially infected nodes to be $o\left(n\right)$, a condition slightly less restrictive than our condition (A3) below. However, none of these previous works have considered multiple layers in the random graph or allowed for activation/deactivation events.

We start by formulating our assumptions in the general case when evolution of the quantities of interest, $\mathit{X}\left(t\right)$, involves a function of the variable $\mathbf{\Theta }\left(t\right)$ defined in the preceding section. In Section 3.2, we state corollaries that relate our result to edge-based models in the special case of static graphs (that is, in the absence of activation or deactivation) [69]. Finally, in Section 3.3 we show that, for a particular class of degree distributions, the evolution of $\mathit{X}\left(t\right)$ decouples from $\mathbf{\Theta }\left(t\right)$. This fact reveals an interesting connection between our limiting system and the one obtained via pair approximation. While others have recently investigated the conditions for exactness of local network moment closures [76,87], we are able to obtain the condition for exactness of a population-level network moment closure. We confirm that, as suggested in [76,86], such a condition depends on the network structure (i.e. degree) heterogeneity. Indeed, in Section 3.3, we employ Theorem 3.1 to prove that the pair approximation approach may or may not give the correct large graph limiting ODEs for the class of LCM stochastic processes described here. In particular, we provide a necessary and sufficient condition on the degree distribution for the two limiting systems to coincide.

3.1. General case

All limits considered below, unless otherwise noted, are with respect to the number of nodes $n\to \mathrm{\infty }$. We use $\stackrel{P}{\to }$ to denote convergence in probability in the product space of the right-continuous with finite left limits (càdlàg) stochastic processes and the space of all random configurations drawn according to an LCM ${\mathcal{G}}_r(\psi ,n)$. We say that a sequence of random variables $Y_n\to \mathrm{\infty }$ with high probability (w.h.p.) if $P(Y_n>k)\to 1$ for any k>0. Let $0<T<\mathrm{\infty }$. We make the following assumptions:

1. For $0<t\le T$, ${\mathit{X}}^{S\text{•}}\left(t\right)\to \mathrm{\infty }$ $w.h.p.$

2. The fractions of initially susceptible, infected, and recovered nodes converge, respectively, to some ${\alpha }^S$, ${\alpha }^I$, ${\alpha }^R$ $\in [0,1]$, i.e. $X^S\left(0\right)/n\stackrel{P}{\to }{\alpha }^S,X^I\left(0\right)/n\stackrel{P}{\to }{\alpha }^I,X^R\left(0\right)/n\stackrel{P}{\to }{\alpha }^R.$ Furthermore, ${\alpha }^S>0$, ${\alpha }^I>0$, and the initially infected and recovered nodes are chosen randomly.

3. $\sum _{\mathit{k}}\parallel \mathit{k}{\parallel }^2{p}_{\mathit{k}}<\mathrm{\infty }.$

The assumption (A1) implies that, for large graphs, some proportion of the individuals remain susceptible on $[0,T]$ and, hence, $\mathbf{\Theta }$ is well-defined in (6(6) $\mathbf{\Theta }\left(t\right)=\exp \left(-\mathit{\beta }{\int }_0^t\frac{{\mathit{X}}^{SI}\left(s\right)}{{\mathit{X}}^{S\text{•}}\left(s\right)}\mathrm{d}s\right)$(6) ). Furthermore, it implies that the average j-degree of a randomly chosen node, i.e. ${\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)$, is positive since $\mathbf{0}<lim inf{n}^{-1}{\mathit{X}}^{S\text{•}}\le \mathbf{\partial }\mathit{\psi }\left(\mathbf{1}\right)$.

Assumption (A2) implies that the initial conditions for the dyads, scaled by n, also converge in probability, i.e. (9) $\begin{array}{rl}& X_j^{SI}\left(0\right)/n\stackrel{P}{\to }{\alpha }_j^{SI},X_j^{\widetilde{SI}}\left(0\right)/n\stackrel{P}{\to }{\alpha }_j^{\widetilde{SI}},X_j^{SS}\left(0\right)/n\stackrel{P}{\to }{\alpha }_j^{SS},\\ & X_j^{\widetilde{SS}}\left(0\right)/n\stackrel{P}{\to }{\alpha }_j^{\widetilde{SS}},j=1,\dots ,r\end{array}$(9) where, for the deactivating layers $j=1,\dots ,m$, ${\alpha }_j^{SI}={\alpha }^S{\alpha }^I{\mu }_j,{\alpha }_j^{\widetilde{SI}}=0,{\alpha }_j^{SS}=({\alpha }^S{)}^2{\mu }_j,{\alpha }_j^{\widetilde{SS}}=0,j=1,\dots ,m$ and, for the activating layers $j=m+1,\dots ,r$, ${\alpha }_j^{SI}=0,{\alpha }_j^{\widetilde{SI}}={\alpha }^S{\alpha }^I{\mu }_j,{\alpha }_j^{SS}=0,{\alpha }_j^{\widetilde{SS}}=({\alpha }^S{)}^2{\mu }_j,j=m+1,\dots ,r.$ To illustrate the above consider, for example, the initial condition ${\alpha }_j^{SI}$ for $j=1,\dots ,m$ as follows. By assumption, all deactivating layer edges are initially activated and all activating layer edges are initially deactivated. Then, according to (A2), the limiting probability of selecting a random node that is susceptible is ${\alpha }^S$. The average number of j-edges a node has is ${\mu }_j$, and the limiting probability that a given edge connects to an infected node is ${\alpha }^I$. Therefore, ${\alpha }_j^{SI}={\alpha }^S{\alpha }^I{\mu }_j$. The other dyad initial conditions are obtained similarly. We denote $\mathit{\alpha }=({\alpha }^S,{\alpha }^I,{\mathit{\alpha }}^{SI},{\mathit{\alpha }}^{\widetilde{SI}},{\mathit{\alpha }}^{SS},{\mathit{\alpha }}^{\widetilde{SS}})$ in what follows.

The assumption (A3) implies that $\sum _{\mathit{k}}{k}_j^2{p}_{\mathit{k}}<\mathrm{\infty }$ for $j=1,\dots ,r$ (i.e. the second moments of the marginal degree distributions are finite) and consequently also that, in each layer, the multigraph constructed by matching half-edges uniformly at random is a simple graph with positive probability (see Section 2 in [43]). Therefore (A3) also guarantees that there is a positive probability of generating a simple LCM graph, i.e. a multilayer graph with a simple graph in each layer.

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 $1\le j,l\le r$, let ${\overline{\kappa }}_{jl}$ be defined by (10) ${\overline{\kappa }}_{jl}\left(\mathit{x}\right)=\frac{\psi \left(\mathit{x}\right){\mathrm{\partial }}_{jl}^2\psi \left(\mathit{x}\right)}{{\mathrm{\partial }}_j\psi \left(\mathit{x}\right){\mathrm{\partial }}_l\psi \left(\mathit{x}\right)}.$(10) As we discuss in Section 3.3.3, ${\overline{\kappa }}_{jl}\left(\mathbf{\Theta }\right)$ can be interpreted as the ratio of the average excess j-degree of a susceptible node chosen randomly as an l-neighbour of an infectious node to the average j-degree of a susceptible node.

We now define the function that in Theorem 3.1 below describes the evolution of $(\mathit{X},\mathbf{\Theta })\left(t\right)$ in the large graph limit. Let $(\mathit{x},\mathit{\theta })=(S,I,\mathbf{[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{I}}\mathbf{]},\mathbf{[}\mathit{S}\mathit{S}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })$ and define $\mathit{H}(\mathit{x},\mathit{\theta })=({\mathit{H}}^{\mathit{X}},{\mathit{H}}^{\mathit{\theta }})(\mathit{x},\mathit{\theta })$ where ${\mathit{H}}^{\mathit{X}}=({\mathcal{H}}^S,{\mathcal{H}}^I,{\mathit{H}}^{SI},{\mathit{H}}^{\widetilde{SI}},{\mathit{H}}^{SS},{\mathit{H}}^{\widetilde{SS}})$ and ${\mathit{H}}^{\mathit{\theta }}$ are given by (11) $\begin{array}{rl}{\mathcal{H}}^S\left(\right[\mathit{S}\mathit{I}\left]\right)& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ {\mathcal{H}}^I(I,[\mathit{S}\mathit{I}\left]\right)& =\sum _{j=1}^r\left({\beta }_r[SI{]}_r\right)-\gamma I,\\ {\mathcal{H}}_j^{SI}\left(\right[\mathit{S}\mathit{I}],[\widetilde{\mathit{S}\mathit{I}}],[\mathit{S}\mathit{S}],\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[SS{]}_j-\left[SI{]}_j\right)\right]\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,\\ {\mathcal{H}}_j^{\widetilde{SI}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{I}}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[\widetilde{SS}{]}_j-\left[\widetilde{SI}{]}_j\right)\right]-(\gamma +{\eta }_j)[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,\\ {\mathcal{H}}_j^{SS}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\mathit{S}\mathit{S}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[SS{]}_j}{S},\\ {\mathcal{H}}_j^{\widetilde{SS}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[\widetilde{SS}{]}_j}{S},\\ {\mathcal{H}}_j^{\theta }\left(\right[\mathit{S}\mathit{I}],\mathit{\theta })& =-{\beta }_j\frac{[SI{]}_j}{{\alpha }^S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)},\end{array}$(11) for $j=1,\dots ,r$.

Theorem 3.1

Strong law of large numbers

Assume conditions (A1)–(A3) for the LCM ${\mathcal{G}}_r(\psi ,n)$. Then, for any $0<T<\mathrm{\infty },$ $\underset{0<t\le T}{sup}\left|\left|(\mathit{X}/n,\mathbf{\Theta })\left(t\right)-(\mathit{x},\mathit{\theta })\left(t\right)\right|\right|\stackrel{P}{\to }0$ where $(\mathit{x},\mathit{\theta })\left(t\right)$ is the solution of (12) $(\mathit{x},\mathit{\theta })\left(t\right)=(\mathit{x},\mathit{\theta })\left(0\right)+{\int }_0^t\mathcal{H}\left(\right(\mathit{x},\mathit{\theta }\left)\right(s\left)\right)\mathrm{d}s$(12) with initial conditions $\mathit{x}\left(0\right)=\mathit{\alpha }$ and $\mathit{\theta }\left(0\right)=\mathbf{1}.$

Theorem 3.1 specifies the large graph limit of the aggregated $SIdaR(r,m)$ process on ${\mathcal{G}}_r(\psi ,n)$ under conditions (A1)–(A3). It says that $(\mathit{X}/n,\mathbf{\Theta }(t\left)\right)$ converges uniformly in probability on any finite interval $[0,T]$ to the solution $(\mathit{x},\mathit{\theta })\left(t\right)$ of the deterministic set of equations given by (12(12) $(\mathit{x},\mathit{\theta })\left(t\right)=(\mathit{x},\mathit{\theta })\left(0\right)+{\int }_0^t\mathcal{H}\left(\right(\mathit{x},\mathit{\theta }\left)\right(s\left)\right)\mathrm{d}s$(12) ). Proof of Theorem 3.1 is provided in Appendix 2. The argument largely follows the standard large volume analysis [50]. The main difficulty is in showing that the terms corresponding to ‘empirical moments’ in system (4(4) $\begin{array}{rlr}{\mathcal{F}}^S\left({\mathit{X}}^{SI}\right)& =-\sum _{l=1}^r{\beta }_l{X}_l^{SI},& \\ {\mathcal{F}}^I(X^I,{\mathit{X}}^{SI})& =\sum _{l=1}^r\left({\beta }_l{X}_l^{SI}\right)-\gamma X^I,& \\ {\mathcal{F}}_j^{SI}({\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SI}},{\mathit{X}}^{SS})& =\sum _{i\in S}\left(\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}(X_j^{SS,i}-X_j^{SI,i})\right)-(\gamma +{\delta }_j)X_j^{SI}+{\eta }_j{X}_j^{\widetilde{SI}},\\ {\mathcal{F}}_j^{\widetilde{SI}}({\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SI}},{\mathit{X}}^{\widetilde{SS}})& =\sum _{i\in S}\left(\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}(X_j^{\widetilde{SS},i}-X_j^{\widetilde{SI},i})\right)-(\gamma +{\eta }_j)X_j^{\widetilde{SI}}+{\delta }_j{X}_j^{SI},\\ {\mathcal{F}}_j^{SS}({\mathit{X}}^{SI},{\mathit{X}}^{SS})& =-2\sum _{i\in S}\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}{X}_j^{SS,i},\\ {\mathcal{F}}_j^{\widetilde{SS}}({\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SS}})& =-2\sum _{i\in S}\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}{X}_j^{\widetilde{SS},i},\end{array}$(4) ) (i.e. those summing over susceptible nodes in the graph and counting exact numbers of neighbours of certain type) can be replaced in the large graph limit (11(11) $\begin{array}{rl}{\mathcal{H}}^S\left(\right[\mathit{S}\mathit{I}\left]\right)& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ {\mathcal{H}}^I(I,[\mathit{S}\mathit{I}\left]\right)& =\sum _{j=1}^r\left({\beta }_r[SI{]}_r\right)-\gamma I,\\ {\mathcal{H}}_j^{SI}\left(\right[\mathit{S}\mathit{I}],[\widetilde{\mathit{S}\mathit{I}}],[\mathit{S}\mathit{S}],\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[SS{]}_j-\left[SI{]}_j\right)\right]\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,\\ {\mathcal{H}}_j^{\widetilde{SI}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{I}}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[\widetilde{SS}{]}_j-\left[\widetilde{SI}{]}_j\right)\right]-(\gamma +{\eta }_j)[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,\\ {\mathcal{H}}_j^{SS}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\mathit{S}\mathit{S}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[SS{]}_j}{S},\\ {\mathcal{H}}_j^{\widetilde{SS}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[\widetilde{SS}{]}_j}{S},\\ {\mathcal{H}}_j^{\theta }\left(\right[\mathit{S}\mathit{I}],\mathit{\theta })& =-{\beta }_j\frac{[SI{]}_j}{{\alpha }^S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)},\end{array}$(11) ) with the corresponding $\overline{\kappa }$ terms encapsulating a population-level average over the heterogeneity in the degrees and excess degrees of susceptible nodes. This may be done due to the properties of the LCM construction with the help of the representation introduced in Remark A.1 in Appendix 2.

3.2. Edge-based 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 ${\delta }_j={\eta }_j=0$ for $j=1,\dots ,n$. Corollary 3.2 below states that in this case our system (12(12) $(\mathit{x},\mathit{\theta })\left(t\right)=(\mathit{x},\mathit{\theta })\left(0\right)+{\int }_0^t\mathcal{H}\left(\right(\mathit{x},\mathit{\theta }\left)\right(s\left)\right)\mathrm{d}s$(12) ) is equivalent to an edge-based model with multiple modes of transmission. This model is the one proposed by Miller and Volz [69] but with a modification to allow for a large number of initially infected nodes (following [67] where Miller modifies the standard SIR edge-based model for such a scenario). In the case that the initially infected nodes are randomly chosen (which we assume in (A2)), the model is given by (13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) The statement of Corollary 3.2 is as follows. Its proof is given in Appendix 3.

Corollary 3.2

Assume ${\delta }_j={\eta }_j=0$ for $j=1,\dots ,n$ and the conditions of Theorem 3.1 hold. Then, the conclusions of Theorem 3.1 hold where $(\mathit{x},\mathit{\theta })\left(t\right)$ is equivalent to the solution of the edge-based model with multiple modes of transmission (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) ).

We further consider the special case of a static, single layer graph (i.e. r=1 and ${\delta }_1={\eta }_1=0$). In the case r=1, (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) ) reduces to the well-known edge-based SIR model of Volz and Miller et al. [70,97], which has been proven to be the large graph limit of the SIR stochastic process on a static configuration model graph [24,43]. By taking r=1 in Corollary 3.2 we have provided an alternative proof of this fact.

Corollary 3.3

Assume $r=1,$ ${\delta }_1={\eta }_1=0$ and the conditions of Theorem 3.1 hold. Then, the conclusions of Theorem 3.1 hold where $(\mathit{x},\theta )\left(t\right)$ is equivalent to the solution of the edge-based SIR model(13(13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) ) with r=1.

3.3. Pairwise limiting systems

In this section we consider a certain class of LCMs where ${\overline{\kappa }}_{ij}$ as defined by Equation (10(10) ${\overline{\kappa }}_{jl}\left(\mathit{x}\right)=\frac{\psi \left(\mathit{x}\right){\mathrm{\partial }}_{jl}^2\psi \left(\mathit{x}\right)}{{\mathrm{\partial }}_j\psi \left(\mathit{x}\right){\mathrm{\partial }}_l\psi \left(\mathit{x}\right)}.$(10) ) is constant. This affords a substantial simplification to the limiting system (11(11) $\begin{array}{rl}{\mathcal{H}}^S\left(\right[\mathit{S}\mathit{I}\left]\right)& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ {\mathcal{H}}^I(I,[\mathit{S}\mathit{I}\left]\right)& =\sum _{j=1}^r\left({\beta }_r[SI{]}_r\right)-\gamma I,\\ {\mathcal{H}}_j^{SI}\left(\right[\mathit{S}\mathit{I}],[\widetilde{\mathit{S}\mathit{I}}],[\mathit{S}\mathit{S}],\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[SS{]}_j-\left[SI{]}_j\right)\right]\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,\\ {\mathcal{H}}_j^{\widetilde{SI}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{I}}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[\widetilde{SS}{]}_j-\left[\widetilde{SI}{]}_j\right)\right]-(\gamma +{\eta }_j)[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,\\ {\mathcal{H}}_j^{SS}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\mathit{S}\mathit{S}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[SS{]}_j}{S},\\ {\mathcal{H}}_j^{\widetilde{SS}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[\widetilde{SS}{]}_j}{S},\\ {\mathcal{H}}_j^{\theta }\left(\right[\mathit{S}\mathit{I}],\mathit{\theta })& =-{\beta }_j\frac{[SI{]}_j}{{\alpha }^S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)},\end{array}$(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. Poisson-type distributions

We define a multivariate Poisson-type (PT) distribution to be a distribution with a pgf that satisfies (14) ${\overline{\kappa }}_{ij}\left(\mathit{x}\right)\equiv {\overline{\kappa }}_{ij}\left(\mathbf{1}\right)={\kappa }_{ij},i,j=1,\dots ,r$(14) where we recall the definition of the normalized average excess degree, ${\kappa }_{ij}$, in Equation (2(2) ${\kappa }_{ji}=\frac{{\mu }_{j|i}^{ex}}{{\mu }_j}=\frac{{\mathrm{\partial }}_{ij}^2\psi \left(\mathbf{1}\right)}{{\mathrm{\partial }}_i\psi \left(\mathbf{1}\right){\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}.$(2) ). At first glance, (14(14) ${\overline{\kappa }}_{ij}\left(\mathit{x}\right)\equiv {\overline{\kappa }}_{ij}\left(\mathbf{1}\right)={\kappa }_{ij},i,j=1,\dots ,r$(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 ${\psi }^{\prime }\left(x\right)={\psi }^{\prime }\left(1\right){\psi }^{\kappa },$ which is satisfied by the univariate Poisson ($\kappa =1$), binomial ($\kappa <1$), and negative binomial ($\kappa >1$) distributions. Note that a k-regular graph (where all nodes have degree k) may be considered a special case of the binomial distribution and that the geometric distribution is a special case of the negative binomial distribution. Bansal et al. [8] have shown that the geometric degree distribution (i.e. the discrete analog of the exponential distribution) gives the best fit for several empirical contact networks. In the multilayer case with r>1, if the marginal degree distributions in different layers are independent, i.e. the pgf $\psi \left(\mathit{x}\right)$ can be written as $\psi \left(\mathit{x}\right)=\prod _{j=1}^r{\psi }_r\left(x_r\right)$, the PT condition (14(14) ${\overline{\kappa }}_{ij}\left(\mathit{x}\right)\equiv {\overline{\kappa }}_{ij}\left(\mathbf{1}\right)={\kappa }_{ij},i,j=1,\dots ,r$(14) ) reduces to the degree distribution in each layer being of (univariate) Poisson-type.

3.3.2. Pairwise model

If the degree distribution is PT as defined by condition (14(14) ${\overline{\kappa }}_{ij}\left(\mathit{x}\right)\equiv {\overline{\kappa }}_{ij}\left(\mathbf{1}\right)={\kappa }_{ij},i,j=1,\dots ,r$(14) ), the limiting system (11(11) $\begin{array}{rl}{\mathcal{H}}^S\left(\right[\mathit{S}\mathit{I}\left]\right)& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ {\mathcal{H}}^I(I,[\mathit{S}\mathit{I}\left]\right)& =\sum _{j=1}^r\left({\beta }_r[SI{]}_r\right)-\gamma I,\\ {\mathcal{H}}_j^{SI}\left(\right[\mathit{S}\mathit{I}],[\widetilde{\mathit{S}\mathit{I}}],[\mathit{S}\mathit{S}],\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[SS{]}_j-\left[SI{]}_j\right)\right]\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,\\ {\mathcal{H}}_j^{\widetilde{SI}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{I}}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[\widetilde{SS}{]}_j-\left[\widetilde{SI}{]}_j\right)\right]-(\gamma +{\eta }_j)[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,\\ {\mathcal{H}}_j^{SS}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\mathit{S}\mathit{S}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[SS{]}_j}{S},\\ {\mathcal{H}}_j^{\widetilde{SS}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[\widetilde{SS}{]}_j}{S},\\ {\mathcal{H}}_j^{\theta }\left(\right[\mathit{S}\mathit{I}],\mathit{\theta })& =-{\beta }_j\frac{[SI{]}_j}{{\alpha }^S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)},\end{array}$(11) ) defining $\mathit{H}=({\mathit{H}}^{\mathit{X}},{\mathit{H}}^{\mathit{\theta }})$ has a particularly simple form. Indeed, substituting the constant ${\kappa }_{jl}$ for ${\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)$ decouples ${\mathit{H}}^{\mathit{X}}$ from ${\mathit{H}}^{\mathit{\theta }}$. We consider the resulting model in this section and introduce some new notation to do so. Let $[XY{]}_j$ and $[\widetilde{XY}{]}_j$ denote, respectively, the number of activated and deactivated edges of type j between a node of status X (either S, I or R) and a node of status Y . Let $[XYZ{]}_{ij}$ denote the number of triples with an active i-edge between nodes of status X and Y and an active j-edge between the node of status Y and one of status Z. Similarly, $[\widetilde{XY}Z{]}_{ij}$ will denote such triples where the i edge is deactivated.

Under the correlation equations approach of Rand [78], triples are needed to describe the evolution of pairs, quadruples (e.g. $[XYZW{]}_{ijk}$) 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 [78]. 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 $p(u=X|A)$ denote the probability that a node u has disease status $X\in \{S,I,R\}$ given an edge (or triple) arrangement A for u. The pair approximation of $[XYZ{]}_{ij}$ is then calculated as follows: (15) $\begin{array}{rl}[XYZ{]}_{ij}& =\left[uYZ{]}_{ij}p\right(u=X|\left[uYZ{]}_{ij}\right)\\ & \approx \left[uYZ{]}_{ij}p\right(u=X|\left[uY{]}_i\right)={\mu }_{i|j}^{ex}[YZ{]}_j\frac{[XY{]}_i}{{\mu }_{i}Y}={\kappa }_{ij}\frac{\left[XY{]}_i\right[YZ{]}_j}{Y}\end{array}$(15) with ${\kappa }_{ij}$ as defined in Equation (2(2) ${\kappa }_{ji}=\frac{{\mu }_{j|i}^{ex}}{{\mu }_j}=\frac{{\mathrm{\partial }}_{ij}^2\psi \left(\mathbf{1}\right)}{{\mathrm{\partial }}_i\psi \left(\mathbf{1}\right){\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}.$(2) ). Note that since alternatively we could have approximated $[XYZ{]}_{ij}\approx [XYu{]}_{ij}p(u=Z|[Yu{]}_j)$ we must have that ${\kappa }_{ji}={\kappa }_{ij}$.

Applying the correlation equations approach using the pair approximation (15(15) $\begin{array}{rl}[XYZ{]}_{ij}& =\left[uYZ{]}_{ij}p\right(u=X|\left[uYZ{]}_{ij}\right)\\ & \approx \left[uYZ{]}_{ij}p\right(u=X|\left[uY{]}_i\right)={\mu }_{i|j}^{ex}[YZ{]}_j\frac{[XY{]}_i}{{\mu }_{i}Y}={\kappa }_{ij}\frac{\left[XY{]}_i\right[YZ{]}_j}{Y}\end{array}$(15) ) to the $SIdaR(r,m)$ dynamics described in Section 2.2 results exactly in the same equations as the limiting system (12(12) $(\mathit{x},\mathit{\theta })\left(t\right)=(\mathit{x},\mathit{\theta })\left(0\right)+{\int }_0^t\mathcal{H}\left(\right(\mathit{x},\mathit{\theta }\left)\right(s\left)\right)\mathrm{d}s$(12) ) in the case of a PT distribution: (16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\sum _{j=1}^r{\beta }_j[SI{]}_j-\gamma I,\\ \frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SI{]}_j\right[SI{]}_l}{S}\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SI}{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SI}{]}_j\right[SI{]}_l}{S}\\ & -({\eta }_j+\gamma )[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[SS{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S},j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SS}{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S},j=1,\dots ,r.\end{array}$(16) Here, we have derived the system (16(16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\sum _{j=1}^r{\beta }_j[SI{]}_j-\gamma I,\\ \frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SI{]}_j\right[SI{]}_l}{S}\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SI}{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SI}{]}_j\right[SI{]}_l}{S}\\ & -({\eta }_j+\gamma )[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[SS{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S},j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SS}{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S},j=1,\dots ,r.\end{array}$(16) ) using absolute pair and triple counts. Notice that, if we scale all variables by the graph size n, the nondimensional quantities satisfy the same system of equations (this holds for any n 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.4. Accordingly, we set the initial conditions to be (17) $(S,I,[SI{]}_1,\dots ,[SI{]}_r,[\widetilde{SI}{]}_1,\dots ,[\widetilde{SI}{]}_r,[SS{]}_1,\dots ,[SS{]}_r,[\widetilde{SS}{]}_1,\dots ,\left[\widetilde{SS}{]}_r\right)\left(0\right)=\mathit{\alpha }$(17) so that they agree with the initial conditions in Theorem 3.1 as $n\to \mathrm{\infty }$.

Observe also that, in fact, we can reduce the dimension of the system (16(16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\sum _{j=1}^r{\beta }_j[SI{]}_j-\gamma I,\\ \frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SI{]}_j\right[SI{]}_l}{S}\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SI}{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SI}{]}_j\right[SI{]}_l}{S}\\ & -({\eta }_j+\gamma )[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[SS{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S},j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SS}{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S},j=1,\dots ,r.\end{array}$(16) ) since we only need to keep track of the deactivated edges for the activating layers, i.e. $j=m+1,\dots ,r$. Also, $[SS{]}_j\equiv 0$ for an activating layer since its initial condition is zero and, hence, we only must track $[SS{]}_j$ for $j=1,\dots ,m$. We refer to system (16(16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\sum _{j=1}^r{\beta }_j[SI{]}_j-\gamma I,\\ \frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SI{]}_j\right[SI{]}_l}{S}\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SI}{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SI}{]}_j\right[SI{]}_l}{S}\\ & -({\eta }_j+\gamma )[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[SS{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S},j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SS}{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S},j=1,\dots ,r.\end{array}$(16) ) with its initial condition (17(17) $(S,I,[SI{]}_1,\dots ,[SI{]}_r,[\widetilde{SI}{]}_1,\dots ,[\widetilde{SI}{]}_r,[SS{]}_1,\dots ,[SS{]}_r,[\widetilde{SS}{]}_1,\dots ,\left[\widetilde{SS}{]}_r\right)\left(0\right)=\mathit{\alpha }$(17) ) as the pairwise model.

The discussion above is summarized in the following corollary.

Corollary 3.4

Assume the conditions of Theorem 3.1 hold for LCM ${\mathcal{G}}_r(\psi ,n)$. Then, the conclusions of Theorem 3.1 hold where $(\mathit{x},\mathit{\theta })\left(t\right)$ is the solution of the pairwise model (16(16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\sum _{j=1}^r{\beta }_j[SI{]}_j-\gamma I,\\ \frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SI{]}_j\right[SI{]}_l}{S}\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SI}{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SI}{]}_j\right[SI{]}_l}{S}\\ & -({\eta }_j+\gamma )[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[SS{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S},j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SS}{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S},j=1,\dots ,r.\end{array}$(16) ) if and only if ${\mathcal{G}}_r(\psi ,n)$ has a multivariate Poisson-type degree distribution.

Furthermore, we can consider the implications of Corollary 3.4 in the static, single layer case. If r=1 and ${\delta }_1={\eta }_1=0$, then the pairwise model (16(16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\sum _{j=1}^r{\beta }_j[SI{]}_j-\gamma I,\\ \frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SI{]}_j\right[SI{]}_l}{S}\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SI}{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SI}{]}_j\right[SI{]}_l}{S}\\ & -({\eta }_j+\gamma )[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[SS{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S},j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SS}{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S},j=1,\dots ,r.\end{array}$(16) ) reduces to the well-known correlation equations model of Rand [78]: (18) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta \left[SI\right],\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta \left[SI\right]-\gamma I,\\ \frac{\mathrm{d}\left[SI\right]}{\mathrm{d}t}& =\beta \kappa \frac{\left[SS\right]\left[SI\right]}{S}-\beta \kappa \frac{[SI{]}^2}{S}-(\beta +\gamma )\left[SI\right],\\ \frac{\mathrm{d}\left[SS\right]}{\mathrm{d}t}& =-2\beta \kappa \frac{\left[SS\right]\left[SI\right]}{S}.\end{array}$(18)

Corollary 3.5

Assume the conditions of Theorem 3.1 hold for ${\mathcal{G}}_1(\psi ,n),$ a static graph $($i.e. ${\delta }_1={\eta }_1=0)$. Then, the conclusions of Theorem 3.1 hold where $(\mathit{x},\theta )\left(t\right)$ is the pairwise SIR model of Rand (18(18) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta \left[SI\right],\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta \left[SI\right]-\gamma I,\\ \frac{\mathrm{d}\left[SI\right]}{\mathrm{d}t}& =\beta \kappa \frac{\left[SS\right]\left[SI\right]}{S}-\beta \kappa \frac{[SI{]}^2}{S}-(\beta +\gamma )\left[SI\right],\\ \frac{\mathrm{d}\left[SS\right]}{\mathrm{d}t}& =-2\beta \kappa \frac{\left[SS\right]\left[SI\right]}{S}.\end{array}$(18) ) if and only if ${\mathcal{G}}_1(\psi ,n)$ has a univariate Poisson-type degree distribution.

Note that together Corollaries 3.2 and 3.4 imply that, in the case of a multivariate PT degree distribution on a static graph, the pairwise model (16(16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\sum _{j=1}^r{\beta }_j[SI{]}_j-\gamma I,\\ \frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SI{]}_j\right[SI{]}_l}{S}\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SI}{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SI}{]}_j\right[SI{]}_l}{S}\\ & -({\eta }_j+\gamma )[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[SS{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S},j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SS}{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S},j=1,\dots ,r.\end{array}$(16) ) is equivalent to the edge-based model with multiple modes of transmission (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) ). Likewise, Corollaries 3.3 and 3.5 indicate that the pairwise SIR model (18(18) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta \left[SI\right],\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta \left[SI\right]-\gamma I,\\ \frac{\mathrm{d}\left[SI\right]}{\mathrm{d}t}& =\beta \kappa \frac{\left[SS\right]\left[SI\right]}{S}-\beta \kappa \frac{[SI{]}^2}{S}-(\beta +\gamma )\left[SI\right],\\ \frac{\mathrm{d}\left[SS\right]}{\mathrm{d}t}& =-2\beta \kappa \frac{\left[SS\right]\left[SI\right]}{S}.\end{array}$(18) ) is equivalent to the edge-based SIR model, (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) ) with r=1, when the distribution is PT. We note that the edge-based SIR model has previously been shown to be equivalent [42,68] to a higher dimensional pairwise model of Eames and Keeling [28]. The latter model stratifies the susceptible nodes by degree and, hence, has dimension K+3 where K represents the number of distinct degrees [68]. The model of dimension K+3 was derived as an approximation to an earlier well-known model of Eames and Keeling, which is of dimension $3K^2/2+3K/2+1$ [28]. We observe that, in fact, (18(18) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\beta \left[SI\right],\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\beta \left[SI\right]-\gamma I,\\ \frac{\mathrm{d}\left[SI\right]}{\mathrm{d}t}& =\beta \kappa \frac{\left[SS\right]\left[SI\right]}{S}-\beta \kappa \frac{[SI{]}^2}{S}-(\beta +\gamma )\left[SI\right],\\ \frac{\mathrm{d}\left[SS\right]}{\mathrm{d}t}& =-2\beta \kappa \frac{\left[SS\right]\left[SI\right]}{S}.\end{array}$(18) ) can be reduced to two differential equations. Separation of variables on $d\left[SS\right]/dS$ (see, e.g. [4] for $\kappa =1$) gives $\left[SS\right]=({\alpha }^S{)}^{2(1-\kappa )}\mu S^{2\kappa }$. Subsequent inspection of $\mathrm{d}\left[SI\right]/\mathrm{d}S$ yields a linear differential equation that can be solved to express $\left[SI\right]$ explicitly as a function of S: (19) $\begin{array}{rl}& =\frac{\beta +\gamma }{\beta (1-\kappa )}S-\mu S^{2\kappa }-\left(\frac{\beta +\gamma }{\beta (1-\kappa )}({\alpha }^S{)}^{1-\kappa }-\mu ({\alpha }^S{)}^{\kappa }-({\alpha }^S{)}^{1-\kappa }{\alpha }^I\mu \right)S^{\kappa },\\ & \kappa \ne 1,\\ \left[SI\right]& =\frac{\beta +\gamma }{\beta }S\log \left(S\right)-\mu S^2-\left(\frac{\beta +\gamma }{\beta }\log \left({\alpha }^S\right)-\mu {\alpha }^S-{\alpha }^I\mu \right)S,\kappa =1.\end{array}$(19) Thus, in Corollary 3.5, we have identified a condition on the degree distribution under which the dimension of a pairwise model that is equivalent to the edge-based SIR model has been reduced from K+3 to two.

3.3.3. Large-graph-consistent pair approximation

Corollary 3.4 and the derivation of the pair approximation (15(15) $\begin{array}{rl}[XYZ{]}_{ij}& =\left[uYZ{]}_{ij}p\right(u=X|\left[uYZ{]}_{ij}\right)\\ & \approx \left[uYZ{]}_{ij}p\right(u=X|\left[uY{]}_i\right)={\mu }_{i|j}^{ex}[YZ{]}_j\frac{[XY{]}_i}{{\mu }_{i}Y}={\kappa }_{ij}\frac{\left[XY{]}_i\right[YZ{]}_j}{Y}\end{array}$(15) ) motivate us to consider a more careful approximation of the triples in the general case when the distribution is not necessarily PT. Let ${\mu }_{i|j}^{ex|SI}$ denote the average excess i-degree of a susceptible node chosen randomly as a j-neighbour of an infectious node and let ${\mu }_i^S$ denote the average i-degree of a susceptible node. We now make more precise our comment in Section 3.1 that ${\overline{\kappa }}_{ij}$ can be interpreted as the limiting ratio of these quantities.

For the LCM ${\mathcal{G}}_r(\psi ,n)$ with $\mathbf{\Theta }$ as defined in Equation (6(6) $\mathbf{\Theta }\left(t\right)=\exp \left(-\mathit{\beta }{\int }_0^t\frac{{\mathit{X}}^{SI}\left(s\right)}{{\mathit{X}}^{S\text{•}}\left(s\right)}\mathrm{d}s\right)$(6) ), we derive ${\mu }_i^S$ and ${\mu }_{i|j}^{ex|SI}$ (see Appendix 4) to be (20) ${\mu }_i^S=\frac{{\mathrm{\Theta }}_i{\mathrm{\partial }}_i\psi \left(\mathbf{\Theta }\right)}{\psi \left(\mathbf{\Theta }\right)}\mathrm{and}{\mu }_{i|j}^{ex|SI}=\frac{{\mathrm{\Theta }}_i{\mathrm{\partial }}_{ij}^2\psi \left(\mathbf{\Theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{\Theta }\right)}.$(20) Thus, we see from the definition of $\overline{\kappa }$ in (10(10) ${\overline{\kappa }}_{jl}\left(\mathit{x}\right)=\frac{\psi \left(\mathit{x}\right){\mathrm{\partial }}_{jl}^2\psi \left(\mathit{x}\right)}{{\mathrm{\partial }}_j\psi \left(\mathit{x}\right){\mathrm{\partial }}_l\psi \left(\mathit{x}\right)}.$(10) ) that ${\overline{\kappa }}_{ij}\left(\mathbf{\Theta }\right)=\frac{{\mu }_{i|j}^{ex|SI}}{{\mu }_i^S}.$ We note that ${\mu }_i^S$ and ${\mu }_{i|j}^{ex|SI}$ are dependent on time and differ, respectively, from ${\mu }_i$ and ${\mu }_{i|j}^{ex}$ defined in Section 2.1. Indeed, susceptible nodes of high degree are preferentially infected [62]. On the other hand, the naive approximation of a triple with ${\kappa }_{ij}={\mu }_{i|j}^{ex}/{\mu }_i$ in (15(15) $\begin{array}{rl}[XYZ{]}_{ij}& =\left[uYZ{]}_{ij}p\right(u=X|\left[uYZ{]}_{ij}\right)\\ & \approx \left[uYZ{]}_{ij}p\right(u=X|\left[uY{]}_i\right)={\mu }_{i|j}^{ex}[YZ{]}_j\frac{[XY{]}_i}{{\mu }_{i}Y}={\kappa }_{ij}\frac{\left[XY{]}_i\right[YZ{]}_j}{Y}\end{array}$(15) ) uses the average degree and excess degree over all nodes, which remain constant for all$t\ge 0$. Note that it is only necessary to approximate triples of the form $[XSI{]}_{ij}$ where $X\in \{S,I\}$. Therefore, we can more carefully derive the pair approximation using ${\mu }_i^S$ and ${\mu }_{i|j}^{ex|SI}$ (21) $\begin{array}{rl}[XSI{]}_{ij}& =\left[uSI{]}_{ij}p\right(u=X|\left[uSI{]}_{ij}\right)\\ & \approx \left[uSI{]}_{ij}p\right(u=X|\left[uS{]}_i\right)={\mu }_{i|j}^{ex|SI}[SI{]}_j\frac{[XS{]}_i}{{\mu }_i^{S}S}={\overline{\kappa }}_{ij}(\mathbf{\Theta })\frac{\left[XS{]}_i\right[SI{]}_j}{S}.\end{array}$(21) Note that the approximation $p(u=X|[uSI{]}_{ij})\approx p(u=X|\left[uS{]}_i\right)$ becomes exact in the large graph limit as the configuration model becomes locally treelike [68,87].

Moreover, it follows from Theorem 3.1 that ${\overline{\kappa }}_{ij}\left(\mathbf{\Theta }\right)\stackrel{P}{\to }{\overline{\kappa }}_{ij}\left(\mathit{\theta }\right)$ uniformly on any finite interval $[0,T]$. We then see from (21(21) $\begin{array}{rl}[XSI{]}_{ij}& =\left[uSI{]}_{ij}p\right(u=X|\left[uSI{]}_{ij}\right)\\ & \approx \left[uSI{]}_{ij}p\right(u=X|\left[uS{]}_i\right)={\mu }_{i|j}^{ex|SI}[SI{]}_j\frac{[XS{]}_i}{{\mu }_i^{S}S}={\overline{\kappa }}_{ij}(\mathbf{\Theta })\frac{\left[XS{]}_i\right[SI{]}_j}{S}.\end{array}$(21) ) that, if we took the correlation equations approach for $SIdaR(r,m)$ dynamics on a general LCM but instead approximated the triples with (22) $[XYZ{]}_{ij}={\overline{\kappa }}_{ij}(\mathit{\theta })\frac{\left[XY{]}_i\right[YZ{]}_j}{Y},$(22) the system of equations obtained would be exactly that given by the limiting system (12(12) $(\mathit{x},\mathit{\theta })\left(t\right)=(\mathit{x},\mathit{\theta })\left(0\right)+{\int }_0^t\mathcal{H}\left(\right(\mathit{x},\mathit{\theta }\left)\right(s\left)\right)\mathrm{d}s$(12) ). In this sense, the pair approximation (22(22) $[XYZ{]}_{ij}={\overline{\kappa }}_{ij}(\mathit{\theta })\frac{\left[XY{]}_i\right[YZ{]}_j}{Y},$(22) ) is ‘correct’, i.e. it is consistent with the large graph limit.

Subsequently, we can view the pairwise system (16(16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\sum _{j=1}^r{\beta }_j[SI{]}_j-\gamma I,\\ \frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SI{]}_j\right[SI{]}_l}{S}\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SI}{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SI}{]}_j\right[SI{]}_l}{S}\\ & -({\eta }_j+\gamma )[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[SS{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S},j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SS}{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S},j=1,\dots ,r.\end{array}$(16) ) as resulting from the approximation (23) ${\overline{\kappa }}_{ij}\left(\mathit{\theta }\right)\approx {\kappa }_{ij}.$(23) For PT distributions, the above relation is exact (‘≈’ may be replace with ‘=’), and the pairwise system is consistent with the large graph limit (i.e. Corollary 3.4). However, for non-PT distributions, the approximation (23(23) ${\overline{\kappa }}_{ij}\left(\mathit{\theta }\right)\approx {\kappa }_{ij}.$(23) ) introduces some error (see Figure 2).

Figure 2. SIR epidemic simulations on static, single-layer graphs with degree distributions given by $\mathrm{Poisson}\left(6\right)$ (blue), a 2-degree distribution in which $p_2=p_9=0.5$ (orange), and a scale-free distribution with $p_k=1.6168k^{-2.01}$ for $1\le k\le 50$ (yellow). (Left) $\overline{\kappa }\left(\theta \right)$ (solid line) is constant and equal to κ (dashed line) for the Poisson distribution but not for the non-PT distributions. (Right) Solutions for the proportion of infected from the large graph limiting system (solid) and the pairwise system (dashed) show discrepancy for the non-PT distributions. Approximate 95% confidence intervals based on 100 stochastic simulations (shaded) and their means (dotted) show agreement with the large graph limiting systems.

We illustrate the results of this section in Figure 2 by simulating an SIR epidemic on three different static, single-layer graphs. The first graph has a Poisson degree distribution with $\mu =6$ (blue). In the second graph, which we refer to as a 2-degree graph, a node has either degree 2 or 9, each with probability 0.5 (orange). The final graph is a so-called scale-free graph [74]; that is, the degree distribution is given by a truncated power law distribution with $p_k=1.6168k^{-2.01}$ for $1\le k\le 50$ (yellow). For each scenario, we simulate the stochastic SIR process on a graph of size $n={10}^5$ and solve the edge-based model (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) ) (which, by Corollary 3.3, is equivalent to the large graph limiting system (12(12) $(\mathit{x},\mathit{\theta })\left(t\right)=(\mathit{x},\mathit{\theta })\left(0\right)+{\int }_0^t\mathcal{H}\left(\right(\mathit{x},\mathit{\theta }\left)\right(s\left)\right)\mathrm{d}s$(12) )) as well as the pairwise system (16(16) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ \frac{\mathrm{d}I}{\mathrm{d}t}& =\sum _{j=1}^r{\beta }_j[SI{]}_j-\gamma I,\\ \frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SI{]}_j\right[SI{]}_l}{S}\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SI}{]}_j}{\mathrm{d}t}& =\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S}-\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SI}{]}_j\right[SI{]}_l}{S}\\ & -({\eta }_j+\gamma )[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,j=1,\dots ,r,\\ \frac{\mathrm{d}[SS{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[SS{]}_j\right[SI{]}_l}{S},j=1,\dots ,r,\\ \frac{\mathrm{d}[\widetilde{SS}{]}_j}{\mathrm{d}t}& =-2\sum _{l=1}^r{\beta }_l{\kappa }_{jl}\frac{\left[\widetilde{SS}{]}_j\right[SI{]}_l}{S},j=1,\dots ,r.\end{array}$(16) ). In Figure 2, we plot the quantity $\overline{\kappa }\left(\theta \right)$ (solid line) as a function of time with the starting value, κ, indicated with a dashed line. The Poisson degree distribution is a PT distribution and, therefore, $\overline{\kappa }\left(\theta \right)\equiv \kappa$; however the latter two graphs do not have PT degree distributions and, indeed, we see a divergence of $\overline{\kappa }\left(\theta \right)$ from κ as the outbreak progresses. This corresponds to a discrepancy between the epidemic curves of the deterministic large graph limiting system (solid) and pairwise system (dashed) in Figure 2. Approximate 95% confidence intervals (shaded) based on 100 stochastic simulations are also shown in Figure 2. For the scale-free graph, $\kappa >\overline{\kappa }\left(\theta \right)$ after the initial time, so using the ‘naive’ approximation (15(15) $\begin{array}{rl}[XYZ{]}_{ij}& =\left[uYZ{]}_{ij}p\right(u=X|\left[uYZ{]}_{ij}\right)\\ & \approx \left[uYZ{]}_{ij}p\right(u=X|\left[uY{]}_i\right)={\mu }_{i|j}^{ex}[YZ{]}_j\frac{[XY{]}_i}{{\mu }_{i}Y}={\kappa }_{ij}\frac{\left[XY{]}_i\right[YZ{]}_j}{Y}\end{array}$(15) ) in the pairwise system leads to overestimating the number of triples. Figure 2 indicates that this consequently results in an overestimate of the incidence of the epidemic (yellow, dashed). In the 2-degree network, on the other hand, κ slightly underestimates $\overline{\kappa }$ and the disease incidence is underestimated by the pairwise system (orange, dashed), although still within the margin of error indicated by the confidence intervals. For the Poisson distribution, the two deterministic curves coincide. The means of the stochastic simulations (dotted) show good agreement with the large graph limiting systems, following Theorem 3.1.

4. Community-healthcare network example

In this section, we discuss the two-layer dynamic network model that was mentioned in the Introduction as a concrete and tractable example of the stochastic SIdaR process. We use this example to illustrate how our Theorem 3.1 can be applied to gain insight into disease dynamics, even in the quite complex setting where network dynamics are tied to infection status. One of the epidemiological issues we aim to understand is how the interplay between edge activation/deactivation and the multilayer network structure affects the ability of disease to invade. In addition, gauging the sensitivity of outbreak size to network parameters and transmissibility along the different edge types would be useful for informing interventions. We demonstrate here that analysis of the pairwise limiting system provides answers to these questions, via application of well-established techniques for analyzing compartmental ODE models of disease dynamics.

Consider the $SIdaR(2,1)$ process, as described in Section 2.2, on a two-layer LCM ${\mathcal{G}}_2(\psi ,n)$ with multivariate PT degree distribution where the two edge types correspond to community and healthcare type contacts. We assume that infected individuals deactivate their community contacts, for example due to decreased mobility or isolation, while they activate their healthcare contacts as they seek care (Figure 1). Note that the healthcare network may include both care provided by healthcare professionals at hospitals or other facilities as well as care provided in the home. This model is motivated by the recent 2014–2015 outbreak of Ebola virus in West Africa. The multitype contact features are particularly relevant for Ebola, given the disproportionate Ebola risk experienced by healthcare workers [22,61] and women (primary caregivers in the home in West Africa) in the 2014 West Africa outbreak as well as Ebola outbreaks in the Democratic Republic of the Congo [49]. There are other aspects of Ebola that are important for transmission, such as an incubation period ranging from 2 to 21 days [36] and disease transmission at funerals [20,36,51]. For simplicity we focus here on community and healthcare transmission and the effect of illness on network structure.

Let C denote community edges and H denote healthcare edges. As in Section 2.2, the stochastic events are assumed to follow independent exponential clocks where the rates of infection along C- and H-edges are, respectively, ${\beta }_C$ and ${\beta }_H$, the rate of deactivation of a C-edge is δ, the rate of activation of an H-edge is η and the rate of recovery is γ.

The $SIdaR(2,1)$ stochastic process above satisfies the conditions of Theorem 3.1 and Corollary 3.4 and, thus, converges to the following system of ODEs in the large graph limit: (24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) As in Section 3.3.2, the variables in (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) are scaled by n and the initial conditions are given by (25) $\begin{array}{rl}S\left(0\right)& ={\alpha }^S,I\left(0\right)={\alpha }^I,R\left(0\right)={\alpha }^R\\ \left[SI{]}_C\right(0)& ={\mu }_C{\alpha }^S{\alpha }^I,\left[SI{]}_H\right(0)=0,[\widetilde{SI}{]}_H\left(0\right)={\mu }_H{\alpha }^S{\alpha }^I\\ \left[SS{]}_C\right(0)& ={\mu }_C({\alpha }^S{)}^2,[SS{]}_H\left(0\right)=0,\left[\widetilde{SS}{]}_H\right(0)={\mu }_H({\alpha }^S{)}^2.\end{array}$(25) Recall that ${\mu }_C$ and ${\mu }_H$ are the average C- and H-degrees, respectively, of a randomly chosen node and that (25(25) $\begin{array}{rl}S\left(0\right)& ={\alpha }^S,I\left(0\right)={\alpha }^I,R\left(0\right)={\alpha }^R\\ \left[SI{]}_C\right(0)& ={\mu }_C{\alpha }^S{\alpha }^I,\left[SI{]}_H\right(0)=0,[\widetilde{SI}{]}_H\left(0\right)={\mu }_H{\alpha }^S{\alpha }^I\\ \left[SS{]}_C\right(0)& ={\mu }_C({\alpha }^S{)}^2,[SS{]}_H\left(0\right)=0,\left[\widetilde{SS}{]}_H\right(0)={\mu }_H({\alpha }^S{)}^2.\end{array}$(25) ) corresponds to all community edges being active and all healthcare edges deactivated at t=0.

Figure 3 compares trajectories of the stochastic $SIdaR(2,1)$ process compared with system (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ). At smaller graph size $n={10}^4$ (Figure 3, left), there are discrepancies between the stochastic trajectories and the limiting ODE (red versus blue curves). As the size of the graph increases, i.e. for $n={10}^5$ (Figure 3, right), the stochastic trajectories are tightly clustered around the ODE solution, illustrating convergence of the stochastic process to the deterministic system given in Theorem 3.1.

Figure 3. Convergence of the stochastic process for infected to its large graph limit ($n\to \mathrm{\infty }$). The degree distribution of the C-layer is $\mathrm{NB}(5,0.706)$ with ${\mu }_C=12$ and ${\kappa }_{CC}=1.2$. The degree distribution of the H-layer is regular, i.e. all nodes have equal degree, with ${\mu }_H=5$ and ${\kappa }_{HH}=0.8$. The size of the graph is given by $n={10}^4$ (left) and $n={10}^5$ (right). The blue shaded regions indicate the approximate 95% confidence intervals based on 100 numerical simulations of the stochastic $SIdaR(2,1)$ processes, conditioned on the large outbreak. The blue lines give the mean of the stochastic simulations, and the red lines show the deterministic solutions to (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ). Initially infected nodes (${\alpha }^I=0.002$, ${\alpha }^S=0.998$) are randomly chosen and other parameters are given by $\gamma =0.08$, $\delta =0.1$, $\eta =0.2$, ${\beta }_C=0.009$, ${\beta }_H=0.015$.

System (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) can now be analyzed to gain insight into how the structure of the different layers of the network and their coupling through activation/deactivation of edges in response to infection affects disease dynamics. In particular, we can compute the basic reproduction number, ${\mathcal{R}}_0$, which determines whether disease invasion is possible [95]. Consider the disease free equilibrium ${\mathit{x}}_{\mathbf{0}}=(1,0,0,0,0,{\mu }_C,{\mu }_H)$. Let $R_C=\frac{{\beta }_C{\kappa }_{CC}{\mu }_C}{{\beta }_C+\gamma +\delta }\mathrm{and}{R}_H=\frac{{\beta }_H{\kappa }_{HH}{\mu }_H\eta }{({\beta }_H+\gamma )(\eta +\gamma )}.$ Note that ${\kappa }_{CC}{\mu }_C={\mu }_C^{ex}$, the average excess C-degree, so we can interpret $R_C$ as the average number of secondary cases transmitted through the community network in a susceptible population and, likewise, $R_H$ represents the secondary cases caused by healthcare transmission. The next-generation matrix method [95] gives (see Appendix A.1) the basic reproduction number as (26a) $\begin{array}{rl}{\mathcal{R}}_0& =\frac{1}{2}\left(R_C+R_H\right)+\frac{1}{2}\sqrt{{\left(R_C+R_H\right)}^2+4R_C{R}_H\left(\frac{{\kappa }_{CH}}{{\kappa }_{CC}}\frac{{\kappa }_{HC}}{{\kappa }_{HH}}-1\right)}\end{array}$(26a) (26b) $\begin{array}{rl}& =\frac{1}{2}\left(R_C+R_H\right)+\frac{1}{2}\sqrt{(R_C-R_H{)}^2+4R_{CH}{R}_{HC}},\end{array}$(26b) where (27) $R_{CH}=\frac{{\beta }_C}{{\beta }_C+\gamma +\delta }{\mu }_{C|H}^{ex}\mathrm{and}{R}_{HC}=\frac{{\beta }_H}{{\beta }_H+\gamma }\frac{\eta }{\eta +\gamma }{\mu }_{H|C}^{ex}.$(27) The $R_{CH}$ term can be interpreted as the number of secondary infections created in the community contact network from an initial infection along an H edge, and $R_{HC}$ as the secondary infections created in the healthcare contact network from an initial infection along a C edge.

Let us comment about the utility of deriving ${\mathcal{R}}_0$ from (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ). The basic reproduction number is a fundamental quantity in epidemiology, and understanding the dependence of ${\mathcal{R}}_0$ on model parameters can be helpful for assessing different intervention strategies. For complex systems such as the stochastic process on a multilayer network considered here, deducing the form for ${\mathcal{R}}_0$ on heuristic grounds can be difficult. Using the next-generation matrix to compute ${\mathcal{R}}_0$ from the limiting system of ODEs takes the guesswork out of this process. As shown in [95], the next generation matrix gives a stability criterion for the disease free equilibrium for the limiting ODEs. Theorem 3.1 shows that the threshold for when the branching process for the stochastic model is subcritical converges to the threshold for the limiting ODEs, i.e. when the ${\mathcal{R}}_0$ expression is equal to one.

Expression (26) is equivalent to the ${\mathcal{R}}_0$ expression found by [7,56], with the difference that the edge activation and deactivation rates enter into the expression. Indeed, the presence of both $R_C$ and $R_H$ in ${\mathcal{R}}_0$ reflects the coupling of the different layers through edge activation/deactivation. In the limit of fast deactivation of community contacts (i.e. $\delta \to \mathrm{\infty }$), ${\mathcal{R}}_0\to R_H$ and disease invasibility depends solely upon the healthcare layer. Similarly, in the limit of slow activation (i.e. $\eta \to 0$), ${\mathcal{R}}_0\to R_C$ and the basic reproduction number is driven by the community layer. For intermediate activation and deactivation rates, ${\mathcal{R}}_0$ depends upon both layers and the multilayer aspect of the model plays an important role in affecting disease dynamics. Both transmission ‘within’ (e.g. $R_C,R_H$ terms) and ‘between’ layers (e.g. $R_{CH},R_{HC}$ terms) contribute to ${\mathcal{R}}_0$, as seen in (26b(26b) $\begin{array}{rl}& =\frac{1}{2}\left(R_C+R_H\right)+\frac{1}{2}\sqrt{(R_C-R_H{)}^2+4R_{CH}{R}_{HC}},\end{array}$(26b) ).

An important point is that (26) shows that knowledge of the reproduction numbers $R_C,R_H$ of the individual layers is not sufficient to determine ${\mathcal{R}}_0$ for the entire network. Instead, ${\mathcal{R}}_0$ depends upon the structure of the different layers through the parameters ${\kappa }_{ij}$. Consider the cross term in (26a(26a) $\begin{array}{rl}{\mathcal{R}}_0& =\frac{1}{2}\left(R_C+R_H\right)+\frac{1}{2}\sqrt{{\left(R_C+R_H\right)}^2+4R_C{R}_H\left(\frac{{\kappa }_{CH}}{{\kappa }_{CC}}\frac{{\kappa }_{HC}}{{\kappa }_{HH}}-1\right)}\end{array}$(26a) ). In the case where ${\kappa }_{CH}{\kappa }_{HC}/{\kappa }_{CC}{\kappa }_{HH}=1$ (for example, independent Poisson degree distributions in each layer), the cross term vanishes and ${\mathcal{R}}_0=R_C+R_H$. In general the sign of the cross term may be either positive or negative, and thus ${\mathcal{R}}_0$ may be either larger or smaller than $R_C+R_H$. In fact for fixed $R_C$ and $R_H$ it is possible for ${\mathcal{R}}_0$ to be either greater than or less than one, depending upon the structure of the layers. This is illustrated in Figure 4 where we plot prevalence (i.e. infected proportion of the population) over the course of an epidemic for two different scenarios for the structure of the healthcare layer. We take the two layer-specific degrees to be independent (${\kappa }_{CH}={\kappa }_{HC}=1$) with the degree distribution of the community layer being Poisson (${\kappa }_{CC}=1$) with ${\mu }_C=10$. We fix $R_C=0.75$ and $R_H=0.5$ and take the degree distribution of the healthcare layer as negative binomial with ${\mu }_H=8$; in the first scenario (blue) we take the healthcare layer degree as $NB(10,4/9)$ (corresponding to ${\kappa }_{HH}=1.1$) while in the second (orange) we take the healthcare layer degree as $NB(1/3,0.96)$ (corresponding to ${\kappa }_{HH}=4$). Figure 4 shows, for each scenario, the deterministic solution to the limiting system (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) (solid line) as well as the approximate 95% empirical confidence interval calculated from 500 stochastic simulations of the corresponding $SIdaR(2,1)$ process with $n={10}^4$ (shaded region). The basic reproduction numbers are calculated from (26a(26a) $\begin{array}{rl}{\mathcal{R}}_0& =\frac{1}{2}\left(R_C+R_H\right)+\frac{1}{2}\sqrt{{\left(R_C+R_H\right)}^2+4R_C{R}_H\left(\frac{{\kappa }_{CH}}{{\kappa }_{CC}}\frac{{\kappa }_{HC}}{{\kappa }_{HH}}-1\right)}\end{array}$(26a) ) to be, respectively, ${\mathcal{R}}_0=1.17$ and ${\mathcal{R}}_0=0.91$. Correspondingly, in the first case a large epidemic occurs while in the second case the initial infection quickly dies out. The increase in ${\kappa }_{HH}$ from the first case to the second corresponds to an increase in the variance of the H degree distribution. In fact, analogous to previous results (see, e.g. [56]), inspection of (26b(26b) $\begin{array}{rl}& =\frac{1}{2}\left(R_C+R_H\right)+\frac{1}{2}\sqrt{(R_C-R_H{)}^2+4R_{CH}{R}_{HC}},\end{array}$(26b) ) reveals that, if ${\mu }_C$ and ${\mu }_H$ are kept constant, ${\mathcal{R}}_0$ is an increasing function of the variances of the C and H degree distributions as well as their covariance.

Figure 4. Prevalence curves for two different scenarios for the structure of a community-healthcare network corresponding to large (top/blue) and small (bottom/orange) outbreaks. The degree distribution of the C-layer is $\mathrm{Pois}\left(10\right)$. The degree distribution of the H-layer is $\mathrm{NB}(10,4/9)$ with ${\kappa }_{HH}=1.1$ in the first case (top/blue) and $\mathrm{NB}(1/3,0.96)$ with ${\kappa }_{HH}=4$ in the second case (bottom/orange). The solid lines show the deterministic solutions to (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) while the shaded regions indicate the approximate 95% confidence intervals based on 500 numerical simulations of the stochastic $SIdaR(2,1)$ processes with $n={10}^4$. Initially infected nodes (${\alpha }^I=0.002$, ${\alpha }^S=0.998$) are randomly chosen and $\gamma =0.1$, $\delta =0.1$, and $\eta =0.3$. We fix $R_C=0.75$ and $R_H=0.5$ which corresponds to ${\beta }_C=0.0162$ and, respectively for the two scenarios, ${\beta }_H=0.0082$ and ${\beta }_H=0.0021$.

In many situations ${\mathcal{R}}_0$ not only determines the ability of disease to invade, but also the size of an outbreak if one occurs [58]. The system of ODEs (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) obtained via Theorem 3.1 and Corollary 3.4 can similarly be analyzed to determine a relation for the final size of the epidemic, as in [4,5]. To illustrate, consider the special case where the degrees in community and healthcare layers are independent with Poisson degree distributions (i.e. ${\kappa }_{CC}={\kappa }_{HH}={\kappa }_{CH}=1$). Let $S_{\mathrm{\infty }}$ denote the fraction of the population that escapes infection. Analysis of a transformed model (as in [4]) or, alternatively, a result of Arino et al. [5] can then be used (see Appendix A.2) to derive the final size relation: (28) $\log \left(\frac{S_{\mathrm{\infty }}}{{\alpha }^S}\right)=-{\mathcal{R}}_0\left({\alpha }^S+{\alpha }^I-S_{\mathrm{\infty }}\right)$(28) which agrees with the classic result for mean-field, homogeneously mixed populations as ${\alpha }^I\to 0$ and ${\alpha }^S\to 1$ [47].

We conclude our analysis of the community-healthcare network model by noting that system (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) is simple enough to be amenable for practical application to outbreaks of interest, for example for parameter estimation and intervention assessment from available data (for example, medical records and contact tracing data). For a given set of parameters system (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) can be used to assess impact of different interventions on preventing outbreaks from occurring (e.g. bringing ${\mathcal{R}}_0<1$) or decreasing the size of an outbreak if it does occur (e.g. by using (28(28) $\log \left(\frac{S_{\mathrm{\infty }}}{{\alpha }^S}\right)=-{\mathcal{R}}_0\left({\alpha }^S+{\alpha }^I-S_{\mathrm{\infty }}\right)$(28) ) to compute the sensitivity of $S_{\mathrm{\infty }}$ to parameters p). Further details on statistical inference and application to specific outbreaks will be presented elsewhere. Here we briefly show how system (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) can be further reduced by finding invariants that allow for dimension reduction. Details for the following methods are provided in Appendix A.3.

Consider the case where ${\kappa }_{CH}\ne {\kappa }_{CC}$. We can then eliminate $[\widetilde{SS}{]}_H$ by expressing it as a function of S and $[SS{]}_C$: (29) $[\widetilde{SS}{]}_H=\frac{{\mu }_H({\alpha }^S{)}^{2(1-\sigma -\lambda )}}{{\mu }_C^{\sigma }}[SS{]}_C^{\sigma }{S}^{2\lambda }$(29) where (30) $\sigma =\frac{{\kappa }_{HH}-{\kappa }_{CH}}{{\kappa }_{CH}-{\kappa }_{CC}}\mathrm{and}\lambda ={\kappa }_{CH}-\sigma {\kappa }_{CC}=\frac{{\kappa }_{CH}^2-{\kappa }_{CC}{\kappa }_{HH}}{{\kappa }_{CH}-{\kappa }_{CC}}.$(30) In the case of independent degrees in different layers (${\kappa }_{CH}=1$) neither of which has a Poisson degree distribution (i.e. ${\kappa }_{CC}\ne 1$ and ${\kappa }_{HH}\ne 1$), we are able to find two additional invariants, Equations (32(32) $\begin{array}{rl}[SI{]}_C& =\frac{{\beta }_C+\gamma +\delta }{{\beta }_C(1-{\kappa }_{CC})}S-[SS{]}_C+\left(\frac{({\beta }_C+\gamma +\delta )}{{\beta }_C({\kappa }_{CC}-1)}+{\mu }_C({\alpha }^S+{\alpha }^I)\right){\left(\frac{[SS{]}_C}{{\mu }_C}\right)}^{1/2}\end{array}$(32) ) and (33(33) $\begin{array}{rl}[\widetilde{SI}{]}_H& =\frac{({\beta }_H+\gamma )(\eta +\gamma )}{\eta {\beta }_H(1-{\kappa }_{HH})}S-\frac{\eta +\gamma }{\eta }[SI{]}_H-[\widetilde{SS}{]}_H\\ & +\left(\frac{({\beta }_H+\gamma )(\eta +\gamma )}{\eta {\beta }_H({\kappa }_{HH}-1)}+{\mu }_H({\alpha }^S+{\alpha }^I)\right){\left(\frac{[\widetilde{SS}{]}_H}{{\mu }_H}\right)}^{1/2}.\end{array}$(33) ) below. The reduced system, of dimension four, is given by (31) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_H[SI{]}_H-{\beta }_C[SI{]}_C\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_H[SI{]}_H+{\beta }_C[SI{]}_C-\gamma I\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H\frac{\left[SS{]}_C\right[SI{]}_H}{S}\end{array}$(31) where $[\widetilde{SS}{]}_H$ is given by (29(29) $[\widetilde{SS}{]}_H=\frac{{\mu }_H({\alpha }^S{)}^{2(1-\sigma -\lambda )}}{{\mu }_C^{\sigma }}[SS{]}_C^{\sigma }{S}^{2\lambda }$(29) ) and (32) $\begin{array}{rl}[SI{]}_C& =\frac{{\beta }_C+\gamma +\delta }{{\beta }_C(1-{\kappa }_{CC})}S-[SS{]}_C+\left(\frac{({\beta }_C+\gamma +\delta )}{{\beta }_C({\kappa }_{CC}-1)}+{\mu }_C({\alpha }^S+{\alpha }^I)\right){\left(\frac{[SS{]}_C}{{\mu }_C}\right)}^{1/2}\end{array}$(32) (33) $\begin{array}{rl}[\widetilde{SI}{]}_H& =\frac{({\beta }_H+\gamma )(\eta +\gamma )}{\eta {\beta }_H(1-{\kappa }_{HH})}S-\frac{\eta +\gamma }{\eta }[SI{]}_H-[\widetilde{SS}{]}_H\\ & +\left(\frac{({\beta }_H+\gamma )(\eta +\gamma )}{\eta {\beta }_H({\kappa }_{HH}-1)}+{\mu }_H({\alpha }^S+{\alpha }^I)\right){\left(\frac{[\widetilde{SS}{]}_H}{{\mu }_H}\right)}^{1/2}.\end{array}$(33) In the case of independent degree layers with Poisson degree distributions (${\kappa }_{CC}={\kappa }_{HH}={\kappa }_{CH}=1$, see Appendix A.2), we have $[SS{]}_C={\mu }_C{S}^2$ and $[\widetilde{SS}{]}_H={\mu }_H{S}^2$. We can further reduce the dimension of the system by one with the following invariant: (34) $\begin{array}{rl}\log \left(\frac{S}{{\alpha }^S}\right)& =-{\mathcal{R}}_0\left({\alpha }^S+{\alpha }^I-S\right)+\frac{{\beta }_C}{{\beta }_C+\gamma +\delta }\frac{[SI{]}_C}{S}\\ & +\frac{{\beta }_H}{{\beta }_H+\gamma }\frac{[SI{]}_H}{S}+\frac{{\beta }_H\eta }{({\beta }_H+\gamma )(\eta +\gamma )}\frac{[\widetilde{SI}{]}_H}{S}.\end{array}$(34)

5. Discussion

The complexity of dynamic multilayer networks makes understanding the disease dynamics evolution on such structures a challenge. Working with simplified models, which are nevertheless capable of retaining most important aspects of network evolution and disease transmission, is essential for gaining biological insight into mechanisms underlying basic disease features such as invasion, persistence, and outbreak size. In this work, we have developed a framework for modeling infectious diseases with multiple modes of transmission in the setting where the network changes in response to infection. Even in this seemingly complicated scenario, it is relatively straightforward to formulate a continuous-time stochastic process by considering transitions in the states of nodes and connected pairs of nodes. However, the state space of the Markov process becomes unmanageable as the size of the network increases and analysis of the non-Markovian aggregate process is complicated. Our main result, Theorem 3.1, rigorously derives the large graph limit of the stochastic process on a layered configuration model graph and, thus, gives a simple model retaining key features of the epidemic process while being amenable to analysis. Moreover, our results extend previous ones for the SIR process on a static, single-layer configuration model network [24,43].

As in [24,43], proof of Theorem 3.1, and expression of the limiting system in the general case, require introduction of the edge-based variable ($\mathbf{\Theta }$). However, in contrast to previous results, we have defined the limiting system here in terms of dyads via the function $\overline{\kappa }$, which also turns out to provide the large-graph-consistent approximation of triples. Consequently, we obtain a simple characterization of the class of degree distributions on random configuration model graphs for which the model obtained via ordinary pair approximation [78] coincides with the limiting system described by Theorem 3.1. The characterization criterion may be formulated in terms of average and excess degrees ($\overline{\kappa }\equiv \kappa$) or, equivalently, in terms of the corresponding identity for the probability generating function. For the single layer case, this condition is satisfied by Poisson but also by binomial and negative binomial distributions. In the non-PT degree distribution case, as our example in Section 3.3.3 illustrates, numerical comparison of κ and $\overline{\kappa }$ could potentially be used to quantify how well a pairwise model accurately reflects the relevant disease dynamics.

Evolving network structure in response to illness is a basic feature that is relevant for many diseases, including Ebola, which was the original motivation for this study. The importance of transmission both within the community and to care-givers for Ebola has empirical support [26,49,80], and these different transmission routes have been incorporated into modeling studies [27,37,51,63,73,79,90,94,101]. Funeral transmission of Ebola is additionally an important concern [41,60], and many modeling studies have incorporated Ebola transmission through unsafe burial practices [51,101,103]. The general multi-layer framework we have presented here with activation and de-activation of edges is flexible and can be extended to include these different Ebola-specific considerations. For example, we can consider a layer corresponding to disease transmission through funerals, with activation of edges in this layer occurring following entry into a death via disease class. Similarly, the basic SIR states for the nodes presented here can be extended to incorporate additional states such an exposed (latent) class, which is important for Ebola [52]. The 2014 West Africa Ebola outbreak prompted an outpouring of modeling studies, utilizing a variety of approaches including ODEs [32,79,90], stochastic processes [27,73], metapopulations [63,94], and contact networks [83,104]. Our work contributes a rigorous approach to understanding how the evolution of network structure in response to infection impacts disease dynamics on layered configuration model networks.

The modeling framework we present is flexible, for example being able accommodate arbitrary joint degree distributions with finite second moment for an arbitrary number of layers. As in the single layer configuration model, however, wiring within each layer occurs at random. Empirical networks often possess features such as community structure [77] or triad closure [100] that are not present in configuration model networks. Despite this, the layered configuration model setting can still yield biological and epidemiological insights due to its tractability for analysis and statistical inference. As demonstrated in the two-layer community-healthcare example in Section 4, the large graph limit we have derived is indeed tractable. The basic reproduction number can be calculated, and its analysis provides insight into how the structure of the different layers and their coupling through edge dynamics affect disease invasion and the final size of the outbreak. Theorem 3.1 and Corollary 3.4 also provide reference points for examining the impact of more realistic network features such as clustering and community structure on disease spread on dynamic multilayer networks.

The general framework presented here can be adapted to diseases with multiple modes of transmission (e.g. those with sexually and non-sexually transmitted infections). Application of this framework to specific diseases, such as Ebola, will require investigation into parameter identifiability and statistical estimation methods. The large graph limit derived here will aid such analysis and, in fact, suggests a hybrid approach in which the node state transitions remain stochastic but the dyads are approximated using the limiting differential equations (or, if possible, invariants such as those found in Section 4). Even for large networks, the resulting Markov process approximation would allow for computationally inexpensive maximum likelihood estimation. Finally, we mention that interventions (e.g. vaccination) can also be incorporated into this framework. The relative simplicity of the limiting system should allow for evaluation of the impact of proposed interventions, for example via sensitivity analysis of ${\mathcal{R}}_0$ and the final outbreak size or via methods of optimal control (e.g. optimize vaccine distribution, see [53]). This could be critical for providing actionable recommendations to public health policymakers with the aim of curbing current epidemics or preventing future outbreaks.

Acknowledgements

We thank Mason Porter and KaYin Leung for their helpful comments during manuscript preparation. We also thank the Mathematical Biosciences Institute at The Ohio State University for its assistance in providing us with space and the necessary computational resources.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Grzegorz A. Rempała http://orcid.org/0000-0002-6307-4555

Additional information

Funding

This research has been supported in part by the Division of Mathematical Sciences at the National Science Foundation under grants DMS-1440386 and RAPID DMS-1513489.

Appendix 1. Summary of notation

The notation for the layered configuration model graph:

n	=	number of nodes
r	=	number of layers in layered configuration model
$\mathit{k}(k_1,\dots ,k_r)$	=	multi-degree with degree $k_i$ in layer i
ψ	=	probability generating function of the multivariate degree distribution
${\mathcal{G}}_r(\psi ,n)$	=	layered configuration model with n nodes, r layers and degree distribution ψ
${\mu }_i$	=	average i-degree
${\kappa }_{ji}$	=	normalized average excess j-degree of an i-neighbour
${\mu }_i^S$	=	average excess i-degree of a susceptible node
${\mu }_{i|j}^{ex|SI}$	=	average excess i-degree of a susceptible node chosen randomly as a j-neighbour of an infectious node

The notation for the $SIdaR(r,m)$ process:

m	=	number of deactivating layers
${\beta }_j$	=	rate of infection along j-edges
${\delta }_j$	=	rate of deactivation (drop) of j-edges
${\eta }_j$	=	rate of activation of j-edges
γ	=	rate of recovery
$X_j^{SI,u}$, $X_j^{SS,u}$	=	number of infectious, susceptible active j-neighbours of u, for susceptible u
$X_j^{\widetilde{SI},u}$, $X_j^{\widetilde{SS},u}$	=	number of infectious, susceptible deactivated j-neighbours of u, for susceptible u
$X^S$, $X^I$, $X^R$	=	number of nodes (of given disease status)
${\mathit{X}}^{SI}$, ${\mathit{X}}^{SS}$	=	number of active edges between nodes (of given status)
${\mathit{X}}^{\widetilde{SI}}$, ${\mathit{X}}^{\widetilde{SS}}$	=	number of deactive edges between nodes (of given status)
${\mathit{X}}^{S\text{•}}$	=	number of edges belonging to susceptible nodes
$\mathbf{\Theta }$	=	quantity defined by Equation (6(6) $\mathbf{\Theta }\left(t\right)=\exp \left(-\mathit{\beta }{\int }_0^t\frac{{\mathit{X}}^{SI}\left(s\right)}{{\mathit{X}}^{S\text{•}}\left(s\right)}\mathrm{d}s\right)$(6) )
${\mathit{F}}^{\mathit{X}}$, ${\mathit{F}}^{\mathbf{\Theta }}$	=	integrable function part of evolution of $\mathit{X}$, $\mathbf{\Theta }$ as in Equations (3(3) $\mathit{X}\left(t\right)=\mathit{X}\left(0\right)+{\int }_0^t{\mathit{F}}^{\mathit{X}}\left(\mathit{X}\right(s\left)\right)\mathrm{d}s+\mathit{M}\left(t\right)$(3) ), (7(7) $\mathbf{\Theta }\left(t\right)=\mathbf{\Theta }\left(0\right)+{\int }_0^t{\mathit{F}}^{\mathbf{\Theta }}\left({\mathit{X}}^{SI}\right(s),{\mathit{X}}^{S\text{•}}(s),\mathbf{\Theta }(s\left)\right)\mathrm{d}s$(7) )

The notation for the large graph limit:

S, I, R	=	number of nodes (of given disease status)
[SI], [SS]	=	number of active edges between nodes (of given status)
[$\widetilde{SI}$], [$\widetilde{SS}$]	=	number of deactive edges between nodes (of given status)
$\mathit{\theta }$	=	large graph limit of $\mathbf{\Theta }$
$\mathit{\alpha }=({\alpha }^S$,…,${\alpha }^{\widetilde{SS}})$	=	initial condition of $\mathit{x}=(S,\dots ,[\widetilde{SS}\left]\right)$, scaled by n
${\mathit{H}}^{\mathit{X}}$, ${\mathit{H}}^{\mathit{\theta }}$	=	integrable function for evolution of $\mathit{X}$, ${\mathit{X}}^{\mathit{\theta }}$ in the large graph limit as in Equation (11(11) $\begin{array}{rl}{\mathcal{H}}^S\left(\right[\mathit{S}\mathit{I}\left]\right)& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ {\mathcal{H}}^I(I,[\mathit{S}\mathit{I}\left]\right)& =\sum _{j=1}^r\left({\beta }_r[SI{]}_r\right)-\gamma I,\\ {\mathcal{H}}_j^{SI}\left(\right[\mathit{S}\mathit{I}],[\widetilde{\mathit{S}\mathit{I}}],[\mathit{S}\mathit{S}],\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[SS{]}_j-\left[SI{]}_j\right)\right]\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,\\ {\mathcal{H}}_j^{\widetilde{SI}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{I}}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[\widetilde{SS}{]}_j-\left[\widetilde{SI}{]}_j\right)\right]-(\gamma +{\eta }_j)[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,\\ {\mathcal{H}}_j^{SS}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\mathit{S}\mathit{S}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[SS{]}_j}{S},\\ {\mathcal{H}}_j^{\widetilde{SS}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[\widetilde{SS}{]}_j}{S},\\ {\mathcal{H}}_j^{\theta }\left(\right[\mathit{S}\mathit{I}],\mathit{\theta })& =-{\beta }_j\frac{[SI{]}_j}{{\alpha }^S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)},\end{array}$(11) )
${\overline{\kappa }}_{jl}$	=	function of network structure used in large graph limit, defined by Equation (10(10) ${\overline{\kappa }}_{jl}\left(\mathit{x}\right)=\frac{\psi \left(\mathit{x}\right){\mathrm{\partial }}_{jl}^2\psi \left(\mathit{x}\right)}{{\mathrm{\partial }}_j\psi \left(\mathit{x}\right){\mathrm{\partial }}_l\psi \left(\mathit{x}\right)}.$(10) )

Appendix 2. Proof of the limit theorem

In this section, we provide the proof of our main result, Theorem 3.1, preceded by two lemmas. The derivation of our results relies on the key observation, summarized in Remark A.1, that in a finite graph the neighbourhood of a susceptible node may be described by a certain multivariate hypergeometric distribution. Lemma A.1 shows that $X^{S_{\mathit{k}}}$ and ${\mathit{X}}^{S\text{•}}$ can be expressed in the limit as functions of $\mathbf{\Theta }$ given by (6(6) $\mathbf{\Theta }\left(t\right)=\exp \left(-\mathit{\beta }{\int }_0^t\frac{{\mathit{X}}^{SI}\left(s\right)}{{\mathit{X}}^{S\text{•}}\left(s\right)}\mathrm{d}s\right)$(6) ). Lemma A.2 shows that the dynamics of the scaled process on the finite graph converges, in the appropriate sense, to the dynamics described by the ODE system (12(12) $(\mathit{x},\mathit{\theta })\left(t\right)=(\mathit{x},\mathit{\theta })\left(0\right)+{\int }_0^t\mathcal{H}\left(\right(\mathit{x},\mathit{\theta }\left)\right(s\left)\right)\mathrm{d}s$(12) ) involving $\mathit{\theta }$. Theorem 3.1 then follows from Lemma A.2 using Doob's and Gronwall's inequalities.

Recall that we take operations on vectors such as multiplication, division, integration and ordering to be coordinatewise. We first provide an important remark about the layer j neighbourhood of a susceptible node of degree $\mathit{k}$, i.e. $\left(X_j^{SI,i}\right(t),X_j^{\widetilde{SI},i}(t),X_j^{SS,i}(t),X_j^{\widetilde{SS},i}(t\left)\right)$ for $i\in S_{\mathit{k}}$. The distribution of the neighbourhood is critical when we consider the expectations of $Q_i^{jl,\left(\mathit{k}\right)}=X_l^{SI,i}(X_j^{SS,i}-X_j^{SI,i})\mathrm{and}{\tilde{Q}}_i^{jl,\left(\mathit{k}\right)}=X_l^{SI,i}(X_j^{\widetilde{SS},i}-X_j^{\widetilde{SI},i}),$ which are mixed moments with respect to the neighbourhood distribution. Recall that we consider the evolution of the $SIdaR(r,m)$ process on a realization of the LCM random graph that has been generated by time t=0. However, we could alternatively consider an equivalent process whereby the graph is revealed dynamically as infections occur, as in [24,43]. In the latter process, a susceptible node $i\in S_{\mathit{k}}$ remains unpaired until it becomes infected. Equivalently we could also pair off all unpaired edges at any time t>0 (uniformly at random and independently in each layer, according to the LCM construction) in order to define the neighbourhood of node i. The following remark is perhaps most easily understood by keeping this equivalent construction in mind and recalling that the probability space considered is that of all random configurations, as described in Section 3.1.

Remark A.1

For $i\in S_{\mathit{k}}$ and $k_j\ge 1$, and conditionally on the process history up to time t, the vector $\left(X_j^{SI,i}\right(t),X_j^{\widetilde{SI},i}(t),X_j^{SS,i}(t),X_j^{\widetilde{SS},i}(t\left)\right)$ follows the multivariate hypergeometric distribution with probability mass function $\begin{array}{rl}P\left(X_j^{SI,i}\right(t)& =n_j^{SI},X_j^{\widetilde{SI},i}\left(t\right)=n_j^{\widetilde{SI}},X_j^{SS,i}\left(t\right)=n_j^{SS},X_j^{\widetilde{SS},i}\left(t\right)=n_j^{\widetilde{SS}})\\ & =\frac{\left(\genfrac{}{}{0em}{}{X_j^{SI}\left(t\right)}{n_j^{SI}}\right)\left(\genfrac{}{}{0em}{}{X_j^{SS}\left(t\right)}{n_j^{SS}}\right)\left(\genfrac{}{}{0em}{}{X_j^{\widetilde{SI}}\left(t\right)}{n_j^{\widetilde{SI}}}\right)\left(\genfrac{}{}{0em}{}{X_j^{\widetilde{SS}}\left(t\right)}{n_j^{\widetilde{SS}}}\right)\left(\genfrac{}{}{0em}{}{X_j^{S\text{•}}\left(t\right)-X_j^{SI}\left(t\right)-X_j^{SS}\left(t\right)-X_j^{\widetilde{SI}}\left(t\right)-X_j^{\widetilde{SS}}\left(t\right)}{k_j-n_j^{SI}-n_j^{SS}-n_j^{\widetilde{SI}}-n_j^{\widetilde{SS}}}\right)}{\left(\genfrac{}{}{0em}{}{X_j^{S\text{•}}\left(t\right)}{k_j}\right)},\end{array}$ supported on the four-simplex $0\le n_j^{SI}+n_j^{\widetilde{SI}}+n_j^{SS}+n_j^{\widetilde{SS}}\le k_j$. This implies $E_h\left[X_j^{SI,i}\right]=k_j{X}_j^{SI}/X_j^{S\text{•}}$ where $E_h$ denotes expectation with respect to the hypergeometric distribution. We also note that, based on the LCM construction, the neighbourhoods of i in distinct layers are independent, i.e. $P\left({\mathit{X}}^{SI,i}\right(t)={\mathit{n}}^{SI},{\mathit{X}}^{SS,i}(t)={\mathit{n}}^{SS})=\prod _{l=1}^{r}P\left(X_l^{SI,i}\right(t)=n_l^{SI},X_l^{SS,i}(t)=n_l^{SS}).$ It follows from the above that the mixed moments are given by (A1) $E_h\left[Q_i^{jl,\left(\mathit{k}\right)}\right]=k_l{k}_j\frac{X_l^{SI}\left(t\right)}{X_l^{S\text{•}}\left(t\right)}\left(\frac{X_j^{SS}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}-\frac{X_j^{SI}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}\right)$(A1) for $l\ne j$ and (A2) $E_h\left[Q_i^{jj,\left(\mathit{k}\right)}\right]=k_j(k_j-1)\frac{X_j^{SI}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}\left(\frac{X_j^{SS}\left(t\right)}{X_j^{S\text{•}}\left(t\right)-1}-\frac{X_j^{SI}\left(t\right)-1}{X_j^{S\text{•}}\left(t\right)-1}\right)-k_j\frac{X_j^{SI}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}.$(A2) Likewise, for any $1\le j,l\le r$, $E_h\left[{\tilde{Q}}_i^{jl,\left(\mathit{k}\right)}\right]=k_l{k}_j\frac{X_l^{SI}\left(t\right)}{X_l^{S\text{•}}\left(t\right)}\left(\frac{X_j^{\widetilde{SS}}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}-\frac{X_j^{\widetilde{SI}}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}\right).$

Keeping in mind this remark, we proceed to prove the first lemma, which shows that $X^{S_{\mathit{k}}}$ and ${\mathit{X}}^{S\text{•}}$ can be expressed in the limit as functions of $\mathbf{\Theta }$.

Lemma A.1

Assume (A1) and $\sum _{\mathit{k}}\parallel \mathit{k}\parallel p_{\mathit{k}}<\mathrm{\infty }$. Then,

1. $\underset{0<t\le T}{sup}|n^{-1}{X}^{S_{\mathit{k}}}-{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}|\stackrel{P}{\to }0$ for any $\mathit{k}\ge 0,$ and

2. $\underset{0<t\le T}{sup}\left|\left|n^{-1}{\mathit{X}}^{S\text{•}}\left(t\right)-{\alpha }^S\mathbf{\Theta }\left(t\right)\mathbf{\partial }\mathit{\psi }\left(\mathbf{\Theta }\right(t\left)\right)\right|\right|\stackrel{P}{\to }0$.

Proof.

$\left(a\right)$. Note $X^{S_{\mathit{k}}}\left(t\right)=\sum _{i=1}^n{Z}_i\left(t\right)$ where $Z_i\left(t\right)\in \{0,1\}$ indicates whether node i is of degree $\mathit{k}$ and susceptible at time t>0. Recall from Remark A.1 that $E{X}_j^{SI,i}=k_j{X}_j^{SI}/X_j^{S\text{•}}$ for $i\in S_{\mathit{k}}$. We claim that (A3) $E{Z}_i\left(t\right)=P(i\in S_{\mathit{k}}(t\left)\right)=n^{-1}S\left(0\right)p_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}.$(A3) Indeed, $E{Z}_i\left(t\right)=P\left(Z_i\right(t)=1)=P(i\in S_{\mathit{k}}(t\left)\right)=P(i\in S_{\mathit{k}}(t)|i\in S_{\mathit{k}}(0\left)\right)P(i\in S_{\mathit{k}}(0\left)\right)$ where $P(i\in S_{\mathit{k}}(0\left)\right)=p_{\mathit{k}}{n}^{-1}S\left(0\right)$ by (A2) and $\begin{array}{rl}P(i\in S_{\mathit{k}}(t)|i\in S_{\mathit{k}}(0\left)\right)& =\exp \left(-\sum _{j=1}^r{\beta }_j{\int }_0^t{E}_h\left[X_j^{SI,i}\left(s\right)\right]\mathrm{d}s\right)\\ & =\exp \left(-\sum _{j=1}^r{\beta }_j{\int }_0^t{k}_j\frac{X_j^{SI}\left(s\right)}{X_j^{S\text{•}}\left(s\right)}\mathrm{d}s\right)={\mathbf{\Theta }}^{\mathit{k}}.\end{array}$ Equation (A3(A3) $E{Z}_i\left(t\right)=P(i\in S_{\mathit{k}}(t\left)\right)=n^{-1}S\left(0\right)p_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}.$(A3) ) then implies that $\left\{X^{S_{\mathit{k}}}\right(t)-S(0\left)p_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\right(t){\}}_{t\ge 0}$ is a càdlàg martingale process with mean zero and finite variation. By the triangle inequality, $\begin{array}{rl}& P\left(\underset{0<t\le T}{sup}\left|n^{-1}{X}^{S_{\mathit{k}}}\left(t\right)-{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)\right|>\epsilon \right)\le P\left(\underset{0<t\le T}{sup}\left|n^{-1}{X}^{S_{\mathit{k}}}\left(t\right)-n^{-1}S\left(0\right)p_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)\right|>\frac{\epsilon }{2}\right)\\ & +P\left(\underset{0<t\le T}{sup}\left|n^{-1}S\left(0\right)p_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)-{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)\right|>\frac{\epsilon }{2}\right).\end{array}$ The second term tends to zero by assumption (A2) and for the first term we have, by Doob's martingale inequality, $P\left(\underset{0<t\le T}{sup}\left|n^{-1}{X}^{S_{\mathit{k}}}\left(t\right)-n^{-1}S\left(0\right)p_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)\right|>\frac{\epsilon }{2}\right)\le {\left(\frac{\epsilon }{2}\right)}^{-2}{n}^{-2}\mathrm{Var}{X}^{S_{\mathit{k}}}\left(T\right).$ Since there are at most n jumps for $X^{S_{\mathit{k}}}$ and each is of size one, it follows that the quadratic variation of $X^{S_{\mathit{k}}}=O\left(n\right)$, which gives the needed result.$\left(b\right)$

(b). By Equation (5(5) ${\mathit{X}}^{S\text{•}}\left(t\right)=\sum _{\mathit{k}}\mathit{k}{X}^{S_{\mathit{k}}}\left(t\right).$(5) ), we can write $n^{-1}{\mathit{X}}^{S\text{•}}\left(t\right)-{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }\left(\mathbf{\Theta }\right(t\left)\right)=\sum _{\mathit{k}}\mathit{k}{n}^{-1}{X}^{S_{\mathit{k}}}\left(t\right)-{\alpha }^S\sum _{\mathit{k}}\mathit{k}{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right).$ Consider arbitrarily large N and $\epsilon >0$. By Markov's inequality, we have $\begin{array}{rl}& P\left(\underset{n}{sup}\underset{0<t\le T}{sup}\left|\left|\sum _{\parallel \mathit{k}\parallel >N}\mathit{k}{n}^{-1}{X}^{S_{\mathit{k}}}\left(t\right)-{\alpha }^S\mathit{k}{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)\right|\right|>\epsilon \right)\\ & \le {\epsilon }^{-1}\sum _{\parallel \mathit{k}\parallel >N}\parallel \mathit{k}\parallel \underset{n}{sup}\underset{0<t\le T}{sup}E\left|n^{-1}{X}^{S_{\mathit{k}}}\left(t\right)-{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)\right|\\ & \le 3{\epsilon }^{-1}\sum _{\parallel \mathit{k}\parallel >N}\parallel \mathit{k}\parallel p_{\mathit{k}}\end{array}$ since $E{X}^{S_{\mathit{k}}}/n\le 2p_{\mathit{k}}$ for n sufficiently large, ${\alpha }^S\le 1$, and we may apply the Monotone Convergence Theorem. So the tail of the sum is negligible since N is arbitrary and, by assumption, $\sum _{\mathit{k}}\parallel \mathit{k}\parallel p_{\mathit{k}}<\mathrm{\infty }$. The result then follows since in $\left(a\right)$ we showed convergence for fixed $\mathit{k}$.

Before proceeding with the next lemma, we give a brief remark on boundedness of our variables and define some useful functions.

Remark A.2

We note that ${\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }\left(\mathbf{\Theta }\right)\le \mathbf{\partial }\mathit{\psi }\left(\mathbf{1}\right)$ and, for sufficiently large n, $n^{-1}{\mathit{X}}^{SI}\le n^{-1}{\mathit{X}}^{S\text{•}}\le 2\mathbf{\partial }\mathit{\psi }\left(\mathbf{1}\right)$ (and likewise for $n^{-1}{\mathit{X}}^{SS}$). By (A1), $n^{-1}{\mathit{X}}^{S\text{•}}$ is bounded away from 0 on $[0,T]$ for finite T and, thus so is $\mathbf{\Theta }$. Furthermore, by Lemma A.1(b), we can take the same lower bound for ${\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }\left(\mathbf{\Theta }\right)$. Let $\xi >0$ be a uniform lower bound for $n^{-1}{\mathit{X}}^{S\text{•}}$, $\mathbf{\Theta }$ and ${\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }\left(\mathbf{\Theta }\right)$. We will use the notation $[\mathit{\xi },2\mathbf{\partial }\mathit{\psi }(\mathbf{1}){]}^r:=[{\xi }_1,2{\mathrm{\partial }}_1\psi \left(\mathbf{1}\right)]\times \cdots \times [{\xi }_r,2{\mathrm{\partial }}_r\psi \left(\mathbf{1}\right)]$ and hence may write ${\mathit{X}}^{S\text{•}}/n\in [\mathit{\xi },2\mathbf{\partial }\mathit{\psi }(\mathbf{1}){]}^r$.

Let ${\mathit{F}}^{\mathit{X}}$ and ${\mathit{F}}^{\mathbf{\Theta }}$ be defined as in (4(4) $\begin{array}{rlr}{\mathcal{F}}^S\left({\mathit{X}}^{SI}\right)& =-\sum _{l=1}^r{\beta }_l{X}_l^{SI},& \\ {\mathcal{F}}^I(X^I,{\mathit{X}}^{SI})& =\sum _{l=1}^r\left({\beta }_l{X}_l^{SI}\right)-\gamma X^I,& \\ {\mathcal{F}}_j^{SI}({\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SI}},{\mathit{X}}^{SS})& =\sum _{i\in S}\left(\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}(X_j^{SS,i}-X_j^{SI,i})\right)-(\gamma +{\delta }_j)X_j^{SI}+{\eta }_j{X}_j^{\widetilde{SI}},\\ {\mathcal{F}}_j^{\widetilde{SI}}({\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SI}},{\mathit{X}}^{\widetilde{SS}})& =\sum _{i\in S}\left(\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}(X_j^{\widetilde{SS},i}-X_j^{\widetilde{SI},i})\right)-(\gamma +{\eta }_j)X_j^{\widetilde{SI}}+{\delta }_j{X}_j^{SI},\\ {\mathcal{F}}_j^{SS}({\mathit{X}}^{SI},{\mathit{X}}^{SS})& =-2\sum _{i\in S}\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}{X}_j^{SS,i},\\ {\mathcal{F}}_j^{\widetilde{SS}}({\mathit{X}}^{SI},{\mathit{X}}^{\widetilde{SS}})& =-2\sum _{i\in S}\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}{X}_j^{\widetilde{SS},i},\end{array}$(4) ) and (8(8) ${\mathit{F}}^{\mathbf{\Theta }}({\mathit{X}}^{SI},{\mathit{X}}^{S\text{•}},\mathbf{\Theta })=-\mathit{\beta }\mathbf{\Theta }{\mathit{X}}^{SI}/{\mathit{X}}^{S\text{•}}.$(8) ) and $\mathit{H}$ as in (11(11) $\begin{array}{rl}{\mathcal{H}}^S\left(\right[\mathit{S}\mathit{I}\left]\right)& =-\sum _{j=1}^r{\beta }_j[SI{]}_j,\\ {\mathcal{H}}^I(I,[\mathit{S}\mathit{I}\left]\right)& =\sum _{j=1}^r\left({\beta }_r[SI{]}_r\right)-\gamma I,\\ {\mathcal{H}}_j^{SI}\left(\right[\mathit{S}\mathit{I}],[\widetilde{\mathit{S}\mathit{I}}],[\mathit{S}\mathit{S}],\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[SS{]}_j-\left[SI{]}_j\right)\right]\\ & -({\beta }_j+\gamma +{\delta }_j)[SI{]}_j+{\eta }_j[\widetilde{SI}{]}_j,\\ {\mathcal{H}}_j^{\widetilde{SI}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{I}}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =\sum _{l=1}^r\left[{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left(\right[\widetilde{SS}{]}_j-\left[\widetilde{SI}{]}_j\right)\right]-(\gamma +{\eta }_j)[\widetilde{SI}{]}_j+{\delta }_j[SI{]}_j,\\ {\mathcal{H}}_j^{SS}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\mathit{S}\mathit{S}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[SS{]}_j}{S},\\ {\mathcal{H}}_j^{\widetilde{SS}}\left(\mathbf{\right[}\mathit{S}\mathit{I}\mathbf{]},\mathbf{[}\widetilde{\mathit{S}\mathit{S}}\mathbf{]},\mathit{\theta })& =-2\sum _{l=1}^r{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{\left[SI{]}_l\right[\widetilde{SS}{]}_j}{S},\\ {\mathcal{H}}_j^{\theta }\left(\right[\mathit{S}\mathit{I}],\mathit{\theta })& =-{\beta }_j\frac{[SI{]}_j}{{\alpha }^S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)},\end{array}$(11) ). Define $\mathbf{\Delta }\left(t\right)=({\mathbf{\Delta }}^{\mathit{X}},{\mathbf{\Delta }}^{\mathbf{\Theta }})\left(t\right)$ where ${\mathbf{\Delta }}^{\mathit{X}}=({\mathrm{\Delta }}^S,{\mathrm{\Delta }}^I,{\mathbf{\Delta }}^{SI},{\mathbf{\Delta }}^{\widetilde{SI}},{\mathbf{\Delta }}^{SS},{\mathbf{\Delta }}^{\widetilde{SS}})$ and ${\mathbf{\Delta }}^{\mathbf{\Theta }}$ are given by (A4) ${\mathbf{\Delta }}^{\mathit{X}}\left(t\right)=n^{-1}{\mathit{F}}^{\mathit{X}}\left(\mathit{X}\right(t\left)\right)-{\mathit{H}}^{\mathit{X}}\left(n^{-1}\mathit{X}\right(t),\mathbf{\Theta }(t\left)\right)$(A4) and (A5) ${\mathbf{\Delta }}^{\mathbf{\Theta }}\left(t\right)={\mathit{F}}^{\mathbf{\Theta }}\left({\mathit{X}}^{SI}\right(t),{\mathit{X}}^{S\text{•}}(t),\mathbf{\Theta }(t\left)\right)-{\mathit{H}}^{\mathit{\theta }}\left(n^{-1}{\mathit{X}}^{SI}\right(t),\mathbf{\Theta }(t\left)\right).$(A5) Lemma A.2 shows that ${\mathbf{\Delta }}^{\mathit{X}}$ and ${\mathbf{\Delta }}^{\mathbf{\Theta }}$ tend to zero uniformly in probability. The convergence of ${\mathbf{\Delta }}^{\mathbf{\Theta }}$ to zero will follow easily from Lemma A.1. The convergence of ${\mathbf{\Delta }}^{\mathit{X}}$ to zero is less obvious and its proof involves consideration of the empirical hypergeometric mixed moments, i.e. $\sum _{i\in S_k}{Q}_i^{jl,\left(\mathit{k}\right)}$ and $\sum _{i\in S_k}{\tilde{Q}}_i^{jl,\left(\mathit{k}\right)}$. We will use the facts that, in the limit, the hypergeometric mixed moments are approximately multinomial and we can replace ${\mathit{X}}^{S\text{•}}$ with a function of $\mathbf{\Theta }$ by Lemma A.1. Therefore, it will be convenient to define the following compensators.

Let $C_h^{jl,\left(\mathit{k}\right)}:[0,T]\to \mathbb{R}$ be given by (A6) $C_h^{jl,\left(\mathit{k}\right)}\left(t\right)=k_l{k}_j\frac{X_l^{SI}\left(t\right)}{X_l^{S\text{•}}\left(t\right)}\left(\frac{X_j^{SS}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}-\frac{X_j^{SI}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}\right),l\ne j,$(A6) and (A7) $C_h^{jj,\left(\mathit{k}\right)}\left(t\right)=k_j(k_j-1)\frac{X_j^{SI}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}\left(\frac{X_j^{SS}\left(t\right)}{X_j^{S\text{•}}\left(t\right)-1}-\frac{X_j^{SI}\left(t\right)-1}{X_j^{S\text{•}}\left(t\right)-1}\right)-k_j\frac{X_j^{SI}\left(t\right)}{X_j^{S\text{•}}\left(t\right)},$(A7) so that the hypergeometric mixed moment in Equation (A1(A1) $E_h\left[Q_i^{jl,\left(\mathit{k}\right)}\right]=k_l{k}_j\frac{X_l^{SI}\left(t\right)}{X_l^{S\text{•}}\left(t\right)}\left(\frac{X_j^{SS}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}-\frac{X_j^{SI}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}\right)$(A1) ) is given by (A8) $E_h\left[Q_i^{jl,\left(\mathit{k}\right)}\right]=C_h^{jl,\left(\mathit{k}\right)}\left(t\right).$(A8) We also define a function related to the multinomial distribution, $C_m^{jl,\left(\mathit{k}\right)}:[0,T]\times [\mathit{\xi },2\mathbf{\partial }\mathit{\psi }(1){]}^r\to \mathbb{R}$, which is given by (A9) $C_m^{jl,\left(\mathit{k}\right)}(t,\mathit{z}(t\left)\right)=k_j{k}_l\frac{n^{-2}{X}_l^{SI}\left(t\right)}{z_l\left(t\right)}\left(\frac{X_j^{SS}\left(t\right)}{z_j\left(t\right)}-\frac{X_j^{SI}\left(t\right)}{z_j\left(t\right)}\right),l\ne j,$(A9) and (A10) $C_m^{jj,\left(\mathit{k}\right)}(t,\mathit{z}(t\left)\right)=k_j(k_j-1)\frac{n^{-2}{X}_j^{SI}}{z_j\left(t\right)}\left(\frac{X_j^{SS}}{z_j\left(t\right)}-\frac{X_j^{SI}}{z_j\left(t\right)}\right)-k_j\frac{n^{-1}{X}_j^{SI}}{z_j\left(t\right)}.$(A10) (Note that if quantities $n^{-1}{X}_j^{SI}/z_j\left(t\right)$ and $n^{-1}{X}_j^{SS}/z_j\left(t\right)$ were actual probabilities then $C_m^{jl,\left(\mathit{k}\right)}(t,\mathit{z}(t\left)\right)$ would be a mixed moment of the multinomial distribution.) Observe that there exists L>0 such that (A11) $C_m^{jl,\left(\mathit{k}\right)}\le L\parallel \mathit{k}{\parallel }^2$(A11) for any j,l (and uniformly in n) sinsce the domain of z is bounded away from 0 and $n^{-1}{\mathit{X}}^{SI}$, $n^{-1}{\mathit{X}}^{SS}$ are uniformly bounded above by Remark A.2. It also follows that $C_m^{jl,\left(\mathit{k}\right)}$ is Lipschitz continuous in $\mathit{z}$.

Lemma A.2

Assume (A1)–(A3). Then,

1. $\underset{0<t\le T}{sup}\parallel {\mathbf{\Delta }}^{\mathbf{\Theta }}\left(t\right)\parallel \stackrel{P}{\to }0,$ and

2. $\underset{0<t\le T}{sup}\parallel {\mathbf{\Delta }}^{\mathit{X}}\left(t\right)\parallel \stackrel{P}{\to }0.$

Proof.

$\left(a\right)$. Define $\mathit{J}:[0,T]\times [\mathit{\xi },2\mathbf{\partial }\mathit{\psi }(\mathbf{1}){]}^r\to {\mathbb{R}}^r$ given by $\mathit{J}(t,\mathit{z}(t\left)\right)=-\beta \mathbf{\Theta }\left(t\right)n^{-1}{\mathit{X}}^{SI}\left(t\right)/\mathit{z}\left(t\right)$. By Remark A.2 and (A1), $(t,n^{-1}{\mathit{X}}^{S\text{•}})$ and $(t,{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }(\mathbf{\Theta }\left)\right)$ are in the domain of $\mathit{J}$ for $t\in [0,T]$. By definition of ${\mathbf{\Delta }}^{\mathbf{\Theta }}$ in (A5(A5) ${\mathbf{\Delta }}^{\mathbf{\Theta }}\left(t\right)={\mathit{F}}^{\mathbf{\Theta }}\left({\mathit{X}}^{SI}\right(t),{\mathit{X}}^{S\text{•}}(t),\mathbf{\Theta }(t\left)\right)-{\mathit{H}}^{\mathit{\theta }}\left(n^{-1}{\mathit{X}}^{SI}\right(t),\mathbf{\Theta }(t\left)\right).$(A5) ), ${\mathbf{\Delta }}^{\mathbf{\Theta }}\left(t\right)=-\beta \frac{{\mathit{X}}^{SI}\left(t\right)\mathbf{\Theta }\left(t\right)}{{\mathit{X}}^{S\text{•}}\left(t\right)}+\beta \frac{n^{-1}{\mathit{X}}^{SI}\left(t\right)}{{\alpha }^S\mathbf{\partial }\mathit{\psi }\left(\mathbf{\Theta }\right(t\left)\right)}=\mathit{J}(t,n^{-1}{\mathit{X}}^{S\text{•}}(t\left)\right)-\mathit{J}(t,{\alpha }^S\mathbf{\Theta }(t\left)\mathbf{\partial }\mathit{\psi }\right(\mathbf{\Theta }\left(t\right)\left)\right).$ Since $\mathit{J}$ is Lipschitz continuous in $\mathit{z}$, $\begin{array}{rl}\underset{0<t\le T}{sup}\parallel {\mathbf{\Delta }}^{\mathbf{\Theta }}\left(t\right)\parallel & =\underset{0<t\le T}{sup}\parallel \mathit{J}(t,n^{-1}{\mathit{X}}^{S\text{•}}(t\left)\right)-\mathit{J}(t,{\alpha }^S\mathbf{\Theta }(t\left)\mathbf{\partial }\mathit{\psi }\right(\mathbf{\Theta }\left(t\right)\left)\right)\parallel \\ & \le L_1\underset{0<t\le T}{sup}\parallel n^{-1}{\mathit{X}}^{S\text{•}}\left(t\right)-{\alpha }^S\mathbf{\Theta }\left(t\right)\mathbf{\partial }\mathit{\psi }\left(\mathbf{\Theta }\right(t\left)\right)\parallel \end{array}$ for some $L_1>0$. The result then follows from Lemma A.1(b) since (A3) implies $\sum _{\mathit{k}}\parallel \mathit{k}\parallel p_{\mathit{k}}<\mathrm{\infty }$.

Regarding part $\left(b\right)$ note that, by definition in (A4(A4) ${\mathbf{\Delta }}^{\mathit{X}}\left(t\right)=n^{-1}{\mathit{F}}^{\mathit{X}}\left(\mathit{X}\right(t\left)\right)-{\mathit{H}}^{\mathit{X}}\left(n^{-1}\mathit{X}\right(t),\mathbf{\Theta }(t\left)\right)$(A4) ), ${\mathrm{\Delta }}^S={\mathrm{\Delta }}^I=0$. We will show that $\underset{0<t\le T}{sup}\parallel {\mathbf{\Delta }}^{SI}\left(t\right)\parallel \stackrel{P}{\to }0$ as it follows similarly for ${\mathbf{\Delta }}^{\widetilde{SI}}$, ${\mathbf{\Delta }}^{SS}$ and ${\mathbf{\Delta }}^{\widetilde{SS}}$ and together these imply $\underset{0<t\le T}{sup}\parallel {\mathbf{\Delta }}^{\mathit{X}}\left(t\right)\parallel \stackrel{P}{\to }0$. In fact, we observe that ${\mathbf{\Delta }}^{SI}=({\mathrm{\Delta }}_1^{SI},\dots ,{\mathrm{\Delta }}_r^{SI})$ and so it suffices to show (A12) $\underset{0<t\le T}{sup}|{\mathrm{\Delta }}_j^{SI}\left(t\right)|\stackrel{P}{\to }0$(A12) for $j=1,\dots ,r$.

Let $1\le j\le r$. We can rewrite ${\mathrm{\Delta }}_j^{SI}$ as $\begin{array}{rl}{\mathrm{\Delta }}_j^{SI}& =n^{-1}{\mathcal{F}}_j^{SI}({\mathit{X}}^{SI},{\mathit{X}}^{SS})-{\mathcal{H}}_j^{SI}(n^{-1}{\mathit{X}}^{SI},n^{-1}{\mathit{X}}^{SS},\mathbf{\Theta })\\ & =n^{-1}\sum _{i\in S}\left(\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}(X_j^{SS,i}-X_j^{SI,i})\right)-\sum _{l=1}^r\left[{\beta }_l{n}^{-2}{X}_l^{SI}(X_j^{SS}-X_j^{SI})\frac{{\mathrm{\partial }}_{jl}^2\psi \left(\mathbf{\Theta }\right)}{{\alpha }^S{\mathrm{\partial }}_j\psi \left(\mathbf{\Theta }\right){\mathrm{\partial }}_l\psi \left(\mathbf{\Theta }\right)}\right]\\ & +{\beta }_j{n}^{-1}{X}_j^{SI}\\ & =\sum _{\mathit{k}}\left[n^{-1}\sum _{i\in S_{\mathit{k}}}\left(\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}(X_j^{SS,i}-X_j^{SI,i})\right)\right.\\ & -\sum _{l\ne j}{\beta }_l{n}^{-2}{X}_l^{SI}\left(X_j^{SS}-X_j^{SI}\right)\frac{k_j{k}_l{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}}{{\alpha }^S{\mathrm{\Theta }}_j{\mathrm{\Theta }}_l{\mathrm{\partial }}_j\psi \left(\mathbf{\Theta }\right){\mathrm{\partial }}_l\psi \left(\mathbf{\Theta }\right)}\\ & \left.+{\beta }_j{n}^{-2}{X}_j^{SI}\left(X_j^{SS}-X_j^{SI}\right)\frac{k_j(k_j-1)p_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}}{{\alpha }^S\left({\mathrm{\Theta }}_j{\mathrm{\partial }}_j\psi \right(\mathbf{\Theta }){)}^2}+{\beta }_j{n}^{-1}{X}_j^{SI}\frac{k_j{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}}{{\mathrm{\Theta }}_j{\mathrm{\partial }}_j\psi \left(\mathbf{\Theta }\right)}\right]\\ & =\sum _{\mathit{k}}\left[n^{-1}\sum _{i\in S_{\mathit{k}}}\left(\sum _{l=1}^r{\beta }_l{X}_l^{SI,i}(X_j^{SS,i}-X_j^{SI,i})\right)-\sum _{l=1}^r{\beta }_l{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}{C}_m^{jl,\left(\mathit{k}\right)}(t,{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }(\mathbf{\Theta }\left)\right)\right].\end{array}$ Thus, we define ${\mathrm{\Delta }}_{jl,\left(\mathit{k}\right)}^{SI}\left(t\right)=n^{-1}\sum _{i\in S_{\mathit{k}}}\left({\beta }_l{X}_l^{SI,i}(X_j^{SS,i}-X_j^{SI,i})\right)-{\beta }_l{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}{C}_m^{jl,\left(\mathit{k}\right)}(t,{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }(\mathbf{\Theta }\left)\right)$ so that ${\mathrm{\Delta }}_j^{SI}=\sum _{l=1}^r\sum _{\mathit{k}}{\mathrm{\Delta }}_{jl,\left(\mathit{k}\right)}^{SI}$. Hence, it suffices to show (A13) $\underset{0<t\le T}{sup}|\sum _{\mathit{k}}{\mathrm{\Delta }}_{jl,\left(\mathit{k}\right)}^{SI}\left(t\right)|\stackrel{P}{\to }0\text{for }j,l=1,\dots ,r.$(A13) This is done in what follows in two separate steps. Consider an arbitrary pair $(j,l)$ $1\le j,\text{ł}\le r.$ We first show that as $N\to \mathrm{\infty }$ (A14) $\underset{n}{sup}\underset{0<t\le T}{sup}\left|\sum _{\parallel \mathit{k}\parallel >N}{\mathrm{\Delta }}_{jl,\left(\mathit{k}\right)}^{SI}\left(t\right)\right|\stackrel{P}{\to }0.$(A14) To this end, observe that (A15) $\begin{array}{rl}& \left|n^{-1}\sum _{\parallel \mathit{k}\parallel >N}\sum _{i\in S_{\mathit{k}}}{\beta }_l{X}_l^{SI,i}\left(t\right)\left(X_j^{SS,i}\right(t)-X_j^{SI,i}(t\left)\right)\right|\\ & \le n^{-1}\sum _{\parallel \mathit{k}\parallel >N}\sum _{i\in S_{\mathit{k}}}|{\beta }_l{k}_l{k}_j|\le n^{-1}{\beta }_l{C}^{\prime }\sum _{\parallel \mathit{k}\parallel >N}\parallel \mathit{k}{\parallel }^2{X}_{S_k}\left(t\right)\le 2{\beta }_l{C}^{\prime }\sum _{\parallel \mathit{k}\parallel >N}\parallel \mathit{k}{\parallel }^2{p}_k\end{array}$(A15) for some $C^{\prime }>0$ since $X_{S_k}/n\le 2p_k$ for n sufficiently large. From the bound on $C_m^{jl,\left(\mathit{k}\right)}$ in (A11(A11) $C_m^{jl,\left(\mathit{k}\right)}\le L\parallel \mathit{k}{\parallel }^2$(A11) ), we also have (A16) $\left|\sum _{\parallel \mathit{k}\parallel >N}{\beta }_l{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}{C}_m^{jl,\left(\mathit{k}\right)}(t,{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }(\mathbf{\Theta }\left)\right)\right|\le {\beta }_{l}L\sum _{\parallel \mathit{k}\parallel >N}\parallel \mathit{k}{\parallel }^2{p}_{\mathit{k}}.$(A16) Then (A14(A14) $\underset{n}{sup}\underset{0<t\le T}{sup}\left|\sum _{\parallel \mathit{k}\parallel >N}{\mathrm{\Delta }}_{jl,\left(\mathit{k}\right)}^{SI}\left(t\right)\right|\stackrel{P}{\to }0.$(A14) ) follows from (A15(A15) $\begin{array}{rl}& \left|n^{-1}\sum _{\parallel \mathit{k}\parallel >N}\sum _{i\in S_{\mathit{k}}}{\beta }_l{X}_l^{SI,i}\left(t\right)\left(X_j^{SS,i}\right(t)-X_j^{SI,i}(t\left)\right)\right|\\ & \le n^{-1}\sum _{\parallel \mathit{k}\parallel >N}\sum _{i\in S_{\mathit{k}}}|{\beta }_l{k}_l{k}_j|\le n^{-1}{\beta }_l{C}^{\prime }\sum _{\parallel \mathit{k}\parallel >N}\parallel \mathit{k}{\parallel }^2{X}_{S_k}\left(t\right)\le 2{\beta }_l{C}^{\prime }\sum _{\parallel \mathit{k}\parallel >N}\parallel \mathit{k}{\parallel }^2{p}_k\end{array}$(A15) ) and (A16(A16) $\left|\sum _{\parallel \mathit{k}\parallel >N}{\beta }_l{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}{C}_m^{jl,\left(\mathit{k}\right)}(t,{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }(\mathbf{\Theta }\left)\right)\right|\le {\beta }_{l}L\sum _{\parallel \mathit{k}\parallel >N}\parallel \mathit{k}{\parallel }^2{p}_{\mathit{k}}.$(A16) ) together with (A3).

Next we show that $\underset{0<t\le T}{sup}|{\mathrm{\Delta }}_{jl,\left(\mathit{k}\right)}^{SI}\left(t\right)|\stackrel{P}{\to }0$ for any $\mathit{k}$. We write (A17) $\begin{array}{rl}|{\mathrm{\Delta }}_{jl,\left(\mathit{k}\right)}^{SI}\left(t\right)|& =\left|n^{-1}\sum _{i\in S_{\mathit{k}}}\left({\beta }_l{X}_l^{SI,i}(X_j^{SS,i}-X_j^{SI,i})\right)-{\beta }_l{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}{C}_m^{jl,\left(\mathit{k}\right)}(t,{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }(\mathbf{\Theta }\left)\right)\right|\\ & \le n^{-1}{\beta }_l\left|\sum _{i\in S_{\mathit{k}}}{Q}_i^{jl,\left(\mathit{k}\right)}-X_{S_k}{C}_h^{jl,\left(\mathit{k}\right)}\left(t\right)\right|\end{array}$(A17) (A18) $+n^{-1}{\beta }_l{X}^{S_{\mathit{k}}}\left|C_h^{jl,\left(\mathit{k}\right)}\left(t\right)-C_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})\right|$(A18) (A19) $+{\beta }_l\left|n^{-1}{X}^{S_{\mathit{k}}}{C}_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})-{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}{C}_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})\right|$(A19) (A20) $+{\beta }_l{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left|C_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})-C_m^{jl,\left(\mathit{k}\right)}(t,{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }(\mathbf{\Theta }\left)\right)\right|$(A20) and we will show that each of these terms tends to zero uniformly in probability.

By Remark A.1 and Equation (A8(A8) $E_h\left[Q_i^{jl,\left(\mathit{k}\right)}\right]=C_h^{jl,\left(\mathit{k}\right)}\left(t\right).$(A8) ), the process $M_h^{jl,\left(\mathit{k}\right)}\left(t\right)=\sum _{i\in S_{\mathit{k}}}{Q}_i^{jl,\left(\mathit{k}\right)}-X_{S_k}{C}_h^{jl,\left(\mathit{k}\right)}\left(t\right)$ is a zero-mean, piecewise-constant càdlàg martingale that jumps only if infection/recovery of a node of degree $\mathit{k}$ (or a neighbour of a node of degree $\mathit{k}$) occurs or activation/drop of a j-edge or l-edge belonging to a node of degree $\mathit{k}$ occurs. Recall that, for each layer, either activation or drops are possible, not both. Consider events impacting a node $u\in S_{\mathit{k}}$. For infection or recovery events, of which there are at most $2(1+k_j+k_l)$ corresponding to infection and recovery of u itself or one of its j- or l-neighbours), the jump size is at most $k_j{k}_l$. For deactivation and activation events of an l or j-edge, of which there are at most $k_l+k_j$ affecting u, the jump size is also at most $k_j{k}_l$. Recall that the number of nodes of degree $\mathit{k}$ is approximately $n{p}_{\mathit{k}}$ for large n. The quadratic variation of $M_h^{jl,\left(\mathit{k}\right)}\left(t\right)$ is the sum of its squared jumps (see, e.g. [4] Chapter 9) and, thus, satisfies $\left[M_h^{jl,\left(\mathit{k}\right)}\right]\left(t\right)=\sum _{s\le t}\left(\delta M_h^{jl,\left(\mathit{k}\right)}\right(s){)}^2\le 2n{p}_{\mathit{k}}(1+k_j+k_l\left)\right(k_j{k}_l{)}^2+n{p}_{\mathit{k}}(k_l+k_j)(k_j{k}_l{)}^2\le L_4\parallel \mathit{k}{\parallel }^{5}n$ for $0<t\le T<\mathrm{\infty }$ and some $L_4>0$. Since $E\left[M_h^{jl,\left(\mathit{k}\right)}\right]\left(t\right)=E\left(M_h^{jl,\left(\mathit{k}\right)}\right(t){)}^2=O(n)$, Doob's martingale inequality implies $\underset{0<t\le T}{sup}\left|n^{-1}{M}_h^{jl,\left(\mathit{k}\right)}\left(t\right)\right|\stackrel{P}{\to }0,$ i.e. the term in (A17(A18) $+n^{-1}{\beta }_l{X}^{S_{\mathit{k}}}\left|C_h^{jl,\left(\mathit{k}\right)}\left(t\right)-C_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})\right|$(A18) ) tends to zero uniformly in probability.

In consideration of the term in (A18(A19) $+{\beta }_l\left|n^{-1}{X}^{S_{\mathit{k}}}{C}_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})-{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}{C}_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})\right|$(A19) ), we note that $n^{-1}{X}_{S_k}\le 1$ and $C_h^{jl,\left(\mathit{k}\right)}\left(t\right)=C_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})$ for $j\ne l$. For the case l=j, we have $\begin{array}{rl}& \underset{0<t\le T}{sup}\left|C_h^{jj,\left(\mathit{k}\right)}\left(t\right)-C_m^{jj,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})\right|=\underset{0<t\le T}{sup}\left|k_j(k_j-1)\frac{X_j^{SI}\left(t\right)}{X_j^{S\text{•}}\left(t\right)}\left(\frac{X_j^{S\text{•}}\left(t\right)+X_j^{SS}\left(t\right)-X_j^{SI}\left(t\right)}{X_j^{S\text{•}}\left(t\right)\left(X_j^{S\text{•}}\right(t)-1)}\right)\right|\\ & \le \frac{2L^{\prime }\parallel \mathit{k}{\parallel }^2}{X_j^{S\text{•}}\left(T\right)-1}\end{array}$ for some $L^{\prime }>0$ and since $X_j^{S\text{•}}\left(t\right)$ is non-increasing on $[0,T]$. Thus, the term in (A18(A18) $+n^{-1}{\beta }_l{X}^{S_{\mathit{k}}}\left|C_h^{jl,\left(\mathit{k}\right)}\left(t\right)-C_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})\right|$(A18) ) tends to zero uniformly in probability by (A1).

For the term in (A19(A19) $+{\beta }_l\left|n^{-1}{X}^{S_{\mathit{k}}}{C}_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})-{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}{C}_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})\right|$(A19) ), we observe $\begin{array}{rl}& \underset{0<t\le T}{sup}\left|n^{-1}{X}^{S_{\mathit{k}}}{C}_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})-{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}{C}_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})\right|\\ & \le L\parallel k{\parallel }^2\underset{0<t\le T}{sup}\left|n^{-1}{X}^{S_{\mathit{k}}}-{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\right|\stackrel{P}{\to }0\end{array}$ by the bound on $C_m^{jl,\left(\mathit{k}\right)}$ in (A11(A11) $C_m^{jl,\left(\mathit{k}\right)}\le L\parallel \mathit{k}{\parallel }^2$(A11) ) and Lemma A.1(a).

Finally, since $C_m^{jl,\left(\mathit{k}\right)}(t,\mathit{z}(t\left)\right)$ is Lipschitz continuous in $\mathit{z}$, we have $\underset{0<t\le T}{sup}\left|C_m^{jl,\left(\mathit{k}\right)}\left(t,n^{-1}{\mathit{X}}^{S\text{•}}\right)-C_m^{jl,\left(\mathit{k}\right)}(t,{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }(\mathbf{\Theta }\left)\right)\right|\le L^{\prime }\underset{0<t\le T}{sup}\left|\left|n^{-1}{\mathit{X}}^{S\text{•}}-{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }\left(\mathbf{\Theta }\right)\right|\right|$ for some $L^{\prime }>0$ and so the term in (A20(A20) $+{\beta }_l{\alpha }^S{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left|C_m^{jl,\left(\mathit{k}\right)}(t,n^{-1}{\mathit{X}}^{S\text{•}})-C_m^{jl,\left(\mathit{k}\right)}(t,{\alpha }^S\mathbf{\Theta }\mathbf{\partial }\mathit{\psi }(\mathbf{\Theta }\left)\right)\right|$(A20) ) tends to zero uniformly in probability by Lemma A.1(b). Therefore, recalling also (A14(A14) $\underset{n}{sup}\underset{0<t\le T}{sup}\left|\sum _{\parallel \mathit{k}\parallel >N}{\mathrm{\Delta }}_{jl,\left(\mathit{k}\right)}^{SI}\left(t\right)\right|\stackrel{P}{\to }0.$(A14) ) we conclude that (A13(A13) $\underset{0<t\le T}{sup}|\sum _{\mathit{k}}{\mathrm{\Delta }}_{jl,\left(\mathit{k}\right)}^{SI}\left(t\right)|\stackrel{P}{\to }0\text{for }j,l=1,\dots ,r.$(A13) ) holds and hence (A12(A12) $\underset{0<t\le T}{sup}|{\mathrm{\Delta }}_j^{SI}\left(t\right)|\stackrel{P}{\to }0$(A12) ) follows.

We may now complete the derivation of Theorem 3.1 via Gronwall's inequality (see, e.g. [4]).

Proof of Theorem 3.1.Remove!

Proof of Theorem 3.1

Recall the definition of $\mathbf{\Delta }=({\mathbf{\Delta }}^{\mathit{X}},{\mathbf{\Delta }}^{\mathbf{\Theta }})$ from Equations (A4(A4) ${\mathbf{\Delta }}^{\mathit{X}}\left(t\right)=n^{-1}{\mathit{F}}^{\mathit{X}}\left(\mathit{X}\right(t\left)\right)-{\mathit{H}}^{\mathit{X}}\left(n^{-1}\mathit{X}\right(t),\mathbf{\Theta }(t\left)\right)$(A4) ) and (A5(A5) ${\mathbf{\Delta }}^{\mathbf{\Theta }}\left(t\right)={\mathit{F}}^{\mathbf{\Theta }}\left({\mathit{X}}^{SI}\right(t),{\mathit{X}}^{S\text{•}}(t),\mathbf{\Theta }(t\left)\right)-{\mathit{H}}^{\mathit{\theta }}\left(n^{-1}{\mathit{X}}^{SI}\right(t),\mathbf{\Theta }(t\left)\right).$(A5) ). Note that, by Equations (3(3) $\mathit{X}\left(t\right)=\mathit{X}\left(0\right)+{\int }_0^t{\mathit{F}}^{\mathit{X}}\left(\mathit{X}\right(s\left)\right)\mathrm{d}s+\mathit{M}\left(t\right)$(3) ) and (7(7) $\mathbf{\Theta }\left(t\right)=\mathbf{\Theta }\left(0\right)+{\int }_0^t{\mathit{F}}^{\mathbf{\Theta }}\left({\mathit{X}}^{SI}\right(s),{\mathit{X}}^{S\text{•}}(s),\mathbf{\Theta }(s\left)\right)\mathrm{d}s$(7) ), $\left(\mathit{X}\right(t)/n,\mathbf{\Theta }(t\left)\right)=\left(\mathit{X}\right(0)/n,\mathbf{\Theta }(0\left)\right)+{\int }_0^t\mathit{H}\left(\mathit{X}\right(s)/n,\mathbf{\Theta }(s\left)\right)\mathrm{d}s+\mathit{E}\left(t\right)$ where $\mathit{E}\left(t\right)=\left(n^{-1}\mathit{M}\right(t),0)+{\int }_0^t\mathbf{\Delta }\left(s\right)\mathrm{d}s,$ where $\mathit{M}\left(t\right)=(M^S,M^I,{\mathit{M}}^{SI},{\mathit{M}}^{\widetilde{SI}},{\mathit{M}}^{SS},{\mathit{M}}^{\widetilde{SS}})\left(t\right)$. Note that each coordinate of $\mathit{M}\left(t\right)$ is a pure jump, càdlàg, zero mean, martingale process. Consider the process $M_j^{SI}$ which, by Equation (3(3) $\mathit{X}\left(t\right)=\mathit{X}\left(0\right)+{\int }_0^t{\mathit{F}}^{\mathit{X}}\left(\mathit{X}\right(s\left)\right)\mathrm{d}s+\mathit{M}\left(t\right)$(3) ), jumps only if infection of a node, recovery of a node or a j-edge drop/activation occurs at time s. Recall that, for each j, either activations or drops are possible, not both. Consider events corresponding to a node of degree $\mathit{k}$. For infection and recovery events the jump size, $\delta M_j^{SI}\left(s\right)$, is not greater than that node's j-degree, and for activation and deactivation events, of which there are at most $2k_j$ affecting that node, the jump size is one. Since the number of nodes of degree $\mathit{k}$ is approximately $n{p}_{\mathit{k}}$ for large n, the corresponding quadratic variation process satisfies $E\left[M_j^{SI}\right]\left(t\right)=E\sum _{s\le t}\left(\delta M_j^{SI}\right(s){)}^2\le 2n\sum _{\mathit{k}}{k}_j^2{p}_{\mathit{k}}+2n\sum _{\mathit{k}}{k}_j{p}_{\mathit{k}}\le 4n\sum _{\mathit{k}}{k}_j^2{p}_{\mathit{k}}=O(n)$ by (A3). Consequently, Doob's martingale inequality implies $\underset{0<t\le T}{sup}\parallel n^{-1}{\mathit{M}}^{SI}\left(t\right)\parallel \stackrel{P}{\to }0$. A similar argument applies also to ${\mathit{M}}^{\widetilde{SI}}$, ${\mathit{M}}^{SS}$, and ${\mathit{M}}^{\widetilde{SS}}$ as well as $M_S$ and $M_I$, both of which make only unit jumps. Since by Lemma A.2 we have $\underset{0<t\le T}{sup}\parallel \mathbf{\Delta }\left(t\right)\parallel \stackrel{P}{\to }0$, this implies $\underset{0<t\le T}{sup}\parallel \mathit{E}\left(t\right)\parallel \stackrel{P}{\to }0$.

Note that $\mathit{H}$ is a (vector valued) Lipschitz continuous function on its domain, which we can take to be $[0,1{]}^2\times ([\mathit{\xi },2\mathbf{\partial }\mathit{\psi }(\mathbf{1}){]}^r{)}^4\times [\mathit{\xi },\mathbf{1}{]}^r$ by Remark A.2. Together with Gronwall's inequality, this implies $\begin{array}{rl}& \underset{0<s\le t}{sup}\parallel (\mathit{X}/n,\mathbf{\Theta })\left(t\right)-(\mathit{x},\mathit{\theta })\left(t\right)\parallel \\ & \le \parallel (\mathit{X}/n,\mathbf{\Theta })\left(0\right)-(\mathit{x},\mathit{\theta })\left(0\right)\parallel +\underset{0<s\le t}{sup}\left|\left|{\int }_0^t\left[\mathit{H}\right((\mathit{X}/n,\mathbf{\Theta })\left(s\right))-\mathit{H}((\mathit{x},\mathit{\theta })\left(s\right)\left)\right]\mathrm{d}s\right|\right|+\underset{0<s\le t}{sup}\parallel \mathit{E}\left(s\right)\parallel \\ & \le \parallel (\mathit{X}/n,\mathbf{\Theta })\left(0\right)-(\mathit{x},\mathit{\theta })\left(0\right)\parallel +L{\int }_0^t\underset{0<s\le t}{sup}\parallel (\mathit{X}/n,\mathbf{\Theta })\left(s\right)-(\mathit{x},\mathit{\theta })\left(s\right)\parallel \mathrm{d}s+\underset{0<s\le t}{sup}\parallel \mathit{E}\left(s\right)\parallel \\ & \le \left(\parallel (\mathit{X}/n,\mathbf{\Theta })\left(0\right)-(\mathit{x},\mathit{\theta })\left(0\right)\parallel +\underset{0<s\le t}{sup}\parallel \mathit{E}\left(s\right)\parallel \right){\mathrm{e}}^{Lt}\stackrel{P}{\to }0,\end{array}$ for some L>0, since the first term in the parenthesis tends to zero by (A2) and (9(9) $\begin{array}{rl}& X_j^{SI}\left(0\right)/n\stackrel{P}{\to }{\alpha }_j^{SI},X_j^{\widetilde{SI}}\left(0\right)/n\stackrel{P}{\to }{\alpha }_j^{\widetilde{SI}},X_j^{SS}\left(0\right)/n\stackrel{P}{\to }{\alpha }_j^{SS},\\ & X_j^{\widetilde{SS}}\left(0\right)/n\stackrel{P}{\to }{\alpha }_j^{\widetilde{SS}},j=1,\dots ,r\end{array}$(9) ). The assertion follows when we take $t=T<\mathrm{\infty }$.

Appendix 3. Equivalence with edge-based multiple modes of transmission model (Proof of Corollary 3.2)

We provide in this section the proof of Corollary 3.2, i.e. equivalence of the system (12(12) $(\mathit{x},\mathit{\theta })\left(t\right)=(\mathit{x},\mathit{\theta })\left(0\right)+{\int }_0^t\mathcal{H}\left(\right(\mathit{x},\mathit{\theta }\left)\right(s\left)\right)\mathrm{d}s$(12) ) with the edge-based model with multiple modes of transmission given by system (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) ). We will use bar notation to denote the variables of the edge-based model.

First, we observe that $S\left(0\right)={\alpha }_S\psi \left(\mathit{\theta }\right(0\left)\right)={\alpha }_S$ and $\frac{\mathrm{d}\left({\alpha }_S\psi \right(\mathit{\theta }\left)\right)}{\mathrm{d}t}={\alpha }_S\sum _j{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=-\sum _j{\beta }_j[SI{]}_j,$ i.e. ${\alpha }_S\psi \left(\mathit{\theta }\right)$ satisfies the same differential equation as S. Therefore, by the uniqueness of the ODE solution, (A21) $S={\alpha }_S\psi \left(\mathit{\theta }\right).$(A21) Secondly, we claim that (A22) $[SS{]}_j={\alpha }_S^2\frac{\left({\mathrm{\partial }}_j\psi \right(\mathit{\theta }){)}^2}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}.$(A22) Indeed, $\left[SS{]}_j\right(0)={\alpha }_S^2{\mathrm{\partial }}_j\psi (\mathbf{1})$ and $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}\left({\alpha }_S^2\frac{\left({\mathrm{\partial }}_j\psi \right(\mathit{\theta }){)}^2}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}\right)& =\frac{2{\alpha }_S^2}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)\\ & =-\frac{2{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}\sum _l{\beta }_l{\mathrm{\partial }}_{jl}^2\psi \left(\mathit{\theta }\right)\frac{[SI{]}_l}{{\mathrm{\partial }}_l\psi \left(\mathit{\theta }\right)}\\ & =-2\sum _l{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\left({\alpha }_S^2\frac{\left({\mathrm{\partial }}_j\psi \right(\mathit{\theta }){)}^2}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}\right),\end{array}$ which is the same differential equation as that satisfied by $[SS{]}_j$, which proves (A22(A22) $[SS{]}_j={\alpha }_S^2\frac{\left({\mathrm{\partial }}_j\psi \right(\mathit{\theta }){)}^2}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}.$(A22) ).

We define ${\phi }^{I,j}={\theta }_j-{\alpha }_S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(1\right)}-\frac{\gamma }{{\beta }_j}(1-{\theta }_j)-{\alpha }_R,j=1,\dots ,r.$ Then, ${\phi }^{I,j}\left(0\right)=1-{\alpha }_S-{\alpha }_R={\alpha }_I$. Moreover, $\begin{array}{rl}\frac{\mathrm{d}{\phi }^{I,j}}{\mathrm{d}t}& =\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}-\frac{{\alpha }_S}{{\mathrm{\partial }}_j\psi \left(1\right)}\sum _l{\mathrm{\partial }}_{jl}^2\psi \left(\mathit{\theta }\right)\frac{\mathrm{d}{\theta }_l}{\mathrm{d}t}+\frac{\gamma }{{\beta }_j}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}\\ & =\frac{1}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}\sum _l{\beta }_l\frac{{\mathrm{\partial }}_{jl}^2\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_l\psi \left(\mathit{\theta }\right)}[SI{]}_l-({\beta }_j+\gamma )\frac{[SI{]}_j}{{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}\\ & =\frac{1}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}\sum _l{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{\psi \left(\mathit{\theta }\right)}[SI{]}_l-({\beta }_j+\gamma )\frac{[SI{]}_j}{{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}.\end{array}$ We now show that $[SI{]}_j/{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)$ satisfies the same differential equation: $\begin{array}{rl}\frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{[SI{]}_j}{{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}\right)& =\frac{\frac{\mathrm{d}[SI{]}_j}{\mathrm{d}t}}{{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}-\frac{\left[SI{]}_j\frac{\mathrm{d}}{\mathrm{d}t}\right({\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right))}{\left({\alpha }_S{\mathrm{\partial }}_j\psi \right(\mathit{\theta }){)}^2}\\ & =\sum _l{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{{\alpha }_{S}S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}\left(\right[SS{]}_j-\left[SI{]}_j\right)-({\beta }_j+\gamma )\frac{[SI{]}_j}{{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}\\ & +\frac{[SI{]}_j}{{\alpha }_S\left({\mathrm{\partial }}_j\psi \right(\mathit{\theta }){)}^2}\sum _l{\beta }_l\frac{{\mathrm{\partial }}_{jl}^2\psi \left(\mathit{\theta }\right)[SI{]}_l}{{\alpha }_S{\mathrm{\partial }}_l\psi \left(\mathit{\theta }\right)}\\ & =\sum _l{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{{\alpha }_{S}S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}\left(\right[SS{]}_j-\left[SI{]}_j\right)-({\beta }_j+\gamma )\frac{[SI{]}_j}{{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}\\ & +\frac{[SI{]}_j}{{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}\sum _l{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{[SI{]}_l}{S}\\ & =\frac{1}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}\sum _l{\beta }_l{\overline{\kappa }}_{jl}\left(\mathit{\theta }\right)\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{\psi \left(\mathit{\theta }\right)}[SI{]}_l-({\beta }_j+\gamma )\frac{[SI{]}_j}{{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}\end{array}$ where we have used (A21(A21) $S={\alpha }_S\psi \left(\mathit{\theta }\right).$(A21) ) and (A22(A22) $[SS{]}_j={\alpha }_S^2\frac{\left({\mathrm{\partial }}_j\psi \right(\mathit{\theta }){)}^2}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}.$(A22) ). Since $\left[SI{]}_j\right(0)/{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right(0\left)\right)={\alpha }_I$ it follows from the uniqueness of the ODE solution that ${\phi }^{I,j}=\frac{[SI{]}_j}{{\alpha }_S{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}.$ Hence, $\mathrm{d}{\theta }_j/\mathrm{d}t=-{\beta }_j{\phi }^{I,j}$, i.e. $\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}=-{\beta }_j{\theta }_j+{\alpha }_S{\beta }_j\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(1\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }_R,$ which is the same differential equation as that for ${\overline{\theta }}_j$ in (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) ). Furthermore, ${\theta }_j\left(1\right)={\overline{\theta }}_j\left(1\right)$ which implies that ${\theta }_j={\overline{\theta }}_j.$ It subsequently follows from (13(13) $\begin{array}{rl}\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}& =-{\beta }_j{\theta }_j+{\beta }_j{\alpha }^S\frac{{\mathrm{\partial }}_j\psi \left(\mathit{\theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{1}\right)}+\gamma (1-{\theta }_j)+{\beta }_j{\alpha }^R,j=1,\dots ,r\\ \frac{\mathrm{d}R}{\mathrm{d}t}& =\gamma I,S={\alpha }^S\psi \left(\mathit{\theta }\right),I=1-S-R,\\ \mathit{\theta }\left(0\right)& =\mathbf{1},R\left(0\right)={\alpha }^R,S\left(0\right)={\alpha }^S,I\left(0\right)=1-{\alpha }^S-{\alpha }^R.\end{array}$(13) ) and Equation (A21(A21) $S={\alpha }_S\psi \left(\mathit{\theta }\right).$(A21) ) that S, I and R are also equivalent for the two models.

Appendix 4. Derivations of ${\mu }_i^S$ and ${\mu }_{i|j}^{ex|SI}$

We provide here the derivations for ${\mu }_i^S$ and ${\mu }_{i|j}^{ex|SI}$ as given in Equation (20(20) ${\mu }_i^S=\frac{{\mathrm{\Theta }}_i{\mathrm{\partial }}_i\psi \left(\mathbf{\Theta }\right)}{\psi \left(\mathbf{\Theta }\right)}\mathrm{and}{\mu }_{i|j}^{ex|SI}=\frac{{\mathrm{\Theta }}_i{\mathrm{\partial }}_{ij}^2\psi \left(\mathbf{\Theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{\Theta }\right)}.$(20) ). Recall that ${\mu }_i^S$ is the average i-degree of a susceptible node. By Equation (A3(A3) $E{Z}_i\left(t\right)=P(i\in S_{\mathit{k}}(t\left)\right)=n^{-1}S\left(0\right)p_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}.$(A3) ), the probability that a node u is susceptible and of degree $\mathit{k}$ is given by $P(u\in S_{\mathit{k}}(t\left)\right)=n^{-1}S\left(0\right)p_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)$. We can then calculate ${\mu }_i^S$ as follows: ${\mu }_i^S\left(t\right)=\sum _{\mathit{k}}{k}_{i}P(u\in S_{\mathit{k}}(t)|u\in S(t\left)\right)=\frac{\sum _{\mathit{k}}{k}_{i}P(u\in S_{\mathit{k}}(t\left)\right)}{\sum _{\mathit{s}}P(u\in S_{\mathit{s}}(t\left)\right)}=\frac{\sum _{\mathit{k}}{k}_i{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)}{\sum _{\mathit{s}}{p}_{\mathit{s}}{\mathbf{\Theta }}^{\mathit{k}}\left(t\right)}=\frac{{\mathrm{\Theta }}_i\left(t\right){\mathrm{\partial }}_i\psi \left(\mathbf{\Theta }\right(t\left)\right)}{\psi \left(\mathbf{\Theta }\right(t\left)\right)}.$ Recall that ${\mu }_{i|j}^{ex|SI}$ is the average excess i-degree of a susceptible node chosen randomly as a j-neighbour of an infectious node. That is, we randomly select a j-edge between a susceptible node and an infectious node, and we calculate the excess i degree of the susceptible node. Let $E_j^{SI}$ denote the set of j-edges between susceptible nodes and infectious nodes and let $E_j^{S_{\mathit{k}}I}$ denote the set of j-edges between susceptible nodes of degree $\mathit{k}$ and infectious nodes, so that $E_j^{SI}={\cup }_{\mathit{k}}{E}_j^{S_{\mathit{k}}I}$. Recall that $X_j^{SI,u}$ is the number of infectious j-neighbours of a susceptible node u. Also recall from Remark A.1 that, given $u\in S_{\mathit{k}}$, the neighbourhood of u has a hypergeometric distribution and $E_h\left[X_j^{SI,u}\right]=k_j{X}_j^{SI}/X_j^{S\text{•}}$. We first assume $i\ne j$ and calculate $\begin{array}{rl}{\mu }_{i|j}^{ex|SI}& =\sum _{\mathit{k}}{k}_{i}P\left(e\in E_j^{S_{\mathit{k}}I}|e\in E_j^{SI}\right)\\ & =\sum _{\mathit{k}}{k}_i\frac{\sum _{l=0}^{k_j}lP\left(u\in S_{\mathit{k}},X_j^{SI,u}=l\right)}{\sum _{\mathit{s}}\sum _{l=0}^{s_j}lP\left(u\in S_{\mathit{s}},X_j^{SI,u}=l\right)}\\ & =\sum _{\mathit{k}}{k}_i\frac{P(u\in S_{\mathit{k}})\sum _{l=0}^{k_j}lP\left(X_j^{SI,u}=l|u\in S_{\mathit{k}}\right)}{\sum _{\mathit{s}}P(u\in S_{\mathit{s}})\sum _{l=0}^{s_j}lP\left(X_j^{SI,u}=l|u\in S_{\mathit{s}}\right)}\\ & =\sum _{\mathit{k}}{k}_i\frac{P(u\in S_{\mathit{k}})k_j{X}_j^{SI}/X_j^{S\text{•}}}{\sum _{\mathit{s}}P(u\in S_{\mathit{s}})s_j{X}_j^{SI}/X_j^{S\text{•}}}\\ & =\frac{\sum _{\mathit{k}}{k}_i{k}_j{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}}{\sum _{\mathit{s}}{s}_j{p}_{\mathit{s}}{\mathbf{\Theta }}^{\mathit{s}}}\\ & =\frac{{\mathrm{\Theta }}_i{\mathrm{\partial }}_{ij}^2\psi \left(\mathbf{\Theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{\Theta }\right)}.\end{array}$ If j=i, we likewise get ${\mu }_{i|i}^{ex|SI}=\frac{\sum _{\mathit{k}}(k_i-1)k_i{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}}{\sum _{\mathit{k}}{k}_i{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}}=\frac{{\mathrm{\Theta }}_i{\mathrm{\partial }}_{ii}^2\psi \left(\mathbf{\Theta }\right)}{{\mathrm{\partial }}_i\psi \left(\mathbf{\Theta }\right)}.$ Therefore, Equation (20(20) ${\mu }_i^S=\frac{{\mathrm{\Theta }}_i{\mathrm{\partial }}_i\psi \left(\mathbf{\Theta }\right)}{\psi \left(\mathbf{\Theta }\right)}\mathrm{and}{\mu }_{i|j}^{ex|SI}=\frac{{\mathrm{\Theta }}_i{\mathrm{\partial }}_{ij}^2\psi \left(\mathbf{\Theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{\Theta }\right)}.$(20) ) follows.

Alternatively, we recall the equivalent model with dynamic graph construction mentioned in Section 2 and consider a j-half edge of an infectious node that is forced to pair with a j-half edge of a susceptible node, which we denote by e, at time t. Let $E_j^S$ denote the set of j-half edges belonging to susceptible nodes and let $E_j^{S_{\mathit{k}}}$ denote those belonging to susceptible nodes of degree $\mathit{k}$. We first assume $i\ne j$ and calculate ${\mu }_{i|j}^{ex|SI}$ as follows: ${\mu }_{i|j}^{ex|SI}=\sum _{\mathit{k}}{k}_{i}P\left(e\in E_j^{S_{\mathit{k}}}|e\in E_j^S\right)=\sum _{\mathit{k}}{k}_i\frac{k_{j}P(u\in S_{\mathit{k}})}{\sum _{\mathit{s}}{s}_{j}P(u\in S_{\mathit{s}})}=\frac{\sum _{\mathit{k}}{k}_i{k}_j{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}}{\sum _{\mathit{k}}{k}_j{p}_{\mathit{k}}{\mathbf{\Theta }}^{\mathit{k}}}=\frac{{\mathrm{\Theta }}_i{\mathrm{\partial }}_{ij}^2\psi \left(\mathbf{\Theta }\right)}{{\mathrm{\partial }}_j\psi \left(\mathbf{\Theta }\right)}.$ The case i=j follows likewise as above.

Appendix 5. Calculations for community-healthcare model

A.1. Basic reproduction number, ${\mathcal{R}}_0$

We use the next-generation matrix method (and corresponding notation) from [95] to calculate ${\mathcal{R}}_0$. We consider I, $[SI{]}_C$, $[SI{]}_H$, and $[\widetilde{SI}{]}_H$ to be the infective compartments (m=4). We consider the disease-free equilibrium for the system (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) given by $x_0=(1,0,0,0,0,{\mu }_C,{\mu }_H)$. The matrix corresponding to terms for new infections, F, is then given by $F=\left[\begin{array}{cccc}0& {\beta }_C& {\beta }_H& 0\\ 0& {\beta }_C{\kappa }_{CC}{\mu }_C& {\beta }_H{\kappa }_{CH}{\mu }_C& 0\\ 0& 0& 0& 0\\ 0& {\beta }_C{\kappa }_{CH}{\mu }_H& {\beta }_H{\kappa }_{HH}{\mu }_H& 0\end{array}\right]$ and the matrix corresponding to terms from all other transitions, V , is $V=\left[\begin{array}{cccc}\gamma & 0& 0& 0\\ 0& {\beta }_C+\gamma +\delta & 0& 0\\ 0& 0& {\beta }_H+\gamma & -\eta \\ 0& 0& 0& \eta +\gamma \end{array}\right].$ Therefore, $V^{-1}=\left[\begin{array}{cccc}1/\gamma & 0& 0& 0\\ 0& 1/({\beta }_C+\gamma +\delta )& 0& 0\\ 0& 0& 1/({\beta }_H+\gamma )& \eta /\left(\right({\beta }_H+\gamma \left)\right(\eta +\gamma \left)\right)\\ 0& 0& 0& 1/(\eta +\gamma )\end{array}\right],$ and the next-generation matrix is given by $\begin{array}{rl}F{V}^{-1}& =\left[\begin{array}{ccc}0& {\beta }_C/({\beta }_C+\gamma +\delta )& {\beta }_H/({\beta }_H+\gamma )\\ 0& {\beta }_C{\kappa }_{CC}{\mu }_C/({\beta }_C+\gamma +\delta )& {\beta }_H{\kappa }_{CH}{\mu }_C/({\beta }_H+\gamma )\\ 0& 0& 0\\ 0& {\beta }_C{\kappa }_{CH}{\mu }_H/({\beta }_C+\gamma +\delta )& {\beta }_H{\kappa }_{HH}{\mu }_H/({\beta }_H+\gamma )\end{array}\right.\\ & \left.\begin{array}{c}{\beta }_H\eta /\left(\right({\beta }_H+\gamma \left)\right(\eta +\gamma \left)\right)\\ {\beta }_H{\kappa }_{CH}{\mu }_C\eta /\left(\right({\beta }_H+\gamma \left)\right(\eta +\gamma \left)\right)\\ 0\\ {\beta }_H{\kappa }_{HH}{\mu }_H\eta /\left(\right({\beta }_H+\gamma \left)\right(\eta +\gamma \left)\right)\end{array}\right].\end{array}$ Then, ${\mathcal{R}}_0$ is the spectral radius of the next-generation matrix, i.e. the largest absolute value of an eigenvalue. We see that $F{V}^{-1}$ has two zero eigenvalues, and the other two eigenvalues are determined by the characteristic polynomial: $p\left(\lambda \right)={\lambda }^2-\left(\frac{{\beta }_C{\kappa }_{CC}{\mu }_C}{{\beta }_C+\gamma +\delta }+\frac{{\beta }_H{\kappa }_{HH}{\mu }_H\eta }{({\beta }_H+\gamma )(\eta +\gamma )}\right)\lambda +\frac{({\kappa }_{CC}{\kappa }_{HH}-{\kappa }_{CH}^2){\beta }_C{\beta }_H{\mu }_C{\mu }_H\eta }{({\beta }_C+\gamma +\delta )({\beta }_H+\gamma )(\eta +\gamma )}.$ Let $R_C=\frac{{\beta }_C{\kappa }_{CC}{\mu }_C}{{\beta }_C+\gamma +\delta },\mathrm{and}{R}_H=\frac{{\beta }_H{\kappa }_{HH}{\mu }_H\eta }{({\beta }_H+\gamma )(\eta +\gamma )}.$ We solve $p\left(\lambda \right)=0$ to determine ${\mathcal{R}}_0=\frac{1}{2}\left(R_C+R_H\right)+\frac{1}{2}\sqrt{{\left(R_C+R_H\right)}^2+4\frac{{\beta }_C{\mu }_C{\beta }_H{\mu }_H\eta }{({\beta }_C+\gamma +\delta )({\beta }_H+\gamma )(\eta +\gamma )}({\kappa }_{CH}^2-{\kappa }_{CC}{\kappa }_{HH})}.$

A.2. Independent Poissons case

We derive an invariant and a final size relation in the case of independent layers with Poisson distributions (i.e. ${\kappa }_{CC}={\kappa }_{HH}={\kappa }_{CH}=1$). First, we observe that, in this case, $\frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}=2\frac{[SS{]}_C}{S}\frac{\mathrm{d}S}{\mathrm{d}t}$ which implies that $[SS{]}_C={\mu }_C{S}^2$. Likewise, $[\widetilde{SS}{]}_H={\mu }_H{S}^2$.

We then transform the system (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) using $Q_i=\frac{[SI{]}_i}{S},{\tilde{Q}}_i=\frac{[\widetilde{SI}{]}_i}{S},\text{for }i=C,H,$ which gives: (A23) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C{Q}_{C}S-{\beta }_H{Q}_{H}S\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C{Q}_{C}S+{\beta }_H{Q}_{H}S-\gamma I\\ \frac{\mathrm{d}{Q}_C}{\mathrm{d}t}& ={\mu }_C({\beta }_C{Q}_{C}S+{\beta }_H{Q}_{H}S)-({\beta }_C+\gamma +\delta )Q_C\\ \frac{\mathrm{d}{Q}_H}{\mathrm{d}t}& =-({\beta }_H+\gamma )Q_H+\eta {\tilde{Q}}_H\\ \frac{\mathrm{d}{\tilde{Q}}_H}{\mathrm{d}t}& ={\mu }_H({\beta }_C{Q}_{C}S+{\beta }_H{Q}_{H}S)-(\eta +\gamma ){\tilde{Q}}_H.\end{array}$(A23) Observe that (A23(A23) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C{Q}_{C}S-{\beta }_H{Q}_{H}S\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C{Q}_{C}S+{\beta }_H{Q}_{H}S-\gamma I\\ \frac{\mathrm{d}{Q}_C}{\mathrm{d}t}& ={\mu }_C({\beta }_C{Q}_{C}S+{\beta }_H{Q}_{H}S)-({\beta }_C+\gamma +\delta )Q_C\\ \frac{\mathrm{d}{Q}_H}{\mathrm{d}t}& =-({\beta }_H+\gamma )Q_H+\eta {\tilde{Q}}_H\\ \frac{\mathrm{d}{\tilde{Q}}_H}{\mathrm{d}t}& ={\mu }_H({\beta }_C{Q}_{C}S+{\beta }_H{Q}_{H}S)-(\eta +\gamma ){\tilde{Q}}_H.\end{array}$(A23) ) fits the form of Arino et al. [5]. Matching their notation, we ignore the decoupled I equation and let the infected compartments be $x=(Q_C,Q_H,{\tilde{Q}}_H{)}^{\mathrm{T}}$, the susceptible compartment y=S, the transmission row vector $\beta b=({\beta }_C,{\beta }_H,0)$, $\mathrm{\Pi }=({\mu }_C,0,{\mu }_H{)}^{\mathrm{T}}$, and $V=\left[\begin{array}{ccc}{\beta }_C+\gamma +\delta & 0& 0\\ 0& {\beta }_H+\gamma & -\eta \\ 0& 0& \eta +\gamma \end{array}\right].$ The system can then be written as (A24) $\frac{\mathrm{d}x}{\mathrm{d}t}=\mathrm{\Pi }y\beta bx-Vx$(A24) (A25) $\frac{\mathrm{d}y}{\mathrm{d}t}=-y\beta bx.$(A25) Recall that $y\left(0\right)=S\left(0\right)={\alpha }_S$, $x\left(0\right)=({\mu }_C{\alpha }_I,0,{\mu }_H{\alpha }_I{)}^{\mathrm{T}}$ and from Equation (26) we have ${\mathcal{R}}_0=\frac{{\beta }_C{\mu }_C}{{\beta }_C+\gamma +\delta }+\frac{{\beta }_H{\mu }_H\eta }{({\beta }_H+\gamma )(\eta +\gamma )}=\beta b{V}^{-1}\mathrm{\Pi }=\beta b{V}^{-1}\frac{x\left(0\right)}{{\alpha }_I}.$ Observe that $\mathrm{\Pi }\frac{\mathrm{d}y}{\mathrm{d}t}+\frac{\mathrm{d}x}{\mathrm{d}t}=-Vx,$ which implies that ${\int }_0^{t}x\left(s\right)\mathrm{d}s=V^{-1}\mathrm{\Pi }\left(y\right(0)-y(t\left)\right)+V^{-1}\left(x\right(0)-x(t\left)\right).$ Then, from integrating Equation (A25(A25) $\frac{\mathrm{d}y}{\mathrm{d}t}=-y\beta bx.$(A25) ), we have $\begin{array}{rl}\log \left(\frac{y\left(t\right)}{y\left(0\right)}\right)& =-\beta b{\int }_0^{t}x\left(s\right)\mathrm{d}s\\ & =-\beta b{V}^{-1}\mathrm{\Pi }\left(y\right(0)-y(t\left)\right)-\beta b{V}^{-1}\left(x\right(0)-x(t\left)\right)\\ & =-{\mathcal{R}}_0({\alpha }_S-y(t\left)\right)-{\mathcal{R}}_0{\alpha }_I+\beta b{V}^{-1}x\left(t\right).\end{array}$ Transforming back to our original variables, this gives the invariant (34(34) $\begin{array}{rl}\log \left(\frac{S}{{\alpha }^S}\right)& =-{\mathcal{R}}_0\left({\alpha }^S+{\alpha }^I-S\right)+\frac{{\beta }_C}{{\beta }_C+\gamma +\delta }\frac{[SI{]}_C}{S}\\ & +\frac{{\beta }_H}{{\beta }_H+\gamma }\frac{[SI{]}_H}{S}+\frac{{\beta }_H\eta }{({\beta }_H+\gamma )(\eta +\gamma )}\frac{[\widetilde{SI}{]}_H}{S}.\end{array}$(34) ): $\begin{array}{rl}\log \left(\frac{S}{{\alpha }_S}\right)& =-{\mathcal{R}}_0\left({\alpha }_S+{\alpha }_I-S\right)+\frac{{\beta }_C}{{\beta }_C+\gamma +\delta }\frac{[SI{]}_C}{S}\\ & +\frac{{\beta }_H}{{\beta }_H+\gamma }\frac{[SI{]}_H}{S}+\frac{{\beta }_H\eta }{({\beta }_H+\gamma )(\eta +\gamma )}\frac{[\widetilde{SI}{]}_H}{S}.\end{array}$ Let $S_{\mathrm{\infty }}$ denote the fraction of the population that escapes infection. Taking the limit $t\to \mathrm{\infty }$ gives the final size relation (28(28) $\log \left(\frac{S_{\mathrm{\infty }}}{{\alpha }^S}\right)=-{\mathcal{R}}_0\left({\alpha }^S+{\alpha }^I-S_{\mathrm{\infty }}\right)$(28) ): $\log \left(\frac{S_{\mathrm{\infty }}}{{\alpha }_S}\right)=-{\mathcal{R}}_0\left({\alpha }_S+{\alpha }_I-S_{\mathrm{\infty }}\right).$

A.3. Model reduction

Let σ and λ be as defined by Equation (30(30) $\sigma =\frac{{\kappa }_{HH}-{\kappa }_{CH}}{{\kappa }_{CH}-{\kappa }_{CC}}\mathrm{and}\lambda ={\kappa }_{CH}-\sigma {\kappa }_{CC}=\frac{{\kappa }_{CH}^2-{\kappa }_{CC}{\kappa }_{HH}}{{\kappa }_{CH}-{\kappa }_{CC}}.$(30) ). We will derive relation (29(29) $[\widetilde{SS}{]}_H=\frac{{\mu }_H({\alpha }^S{)}^{2(1-\sigma -\lambda )}}{{\mu }_C^{\sigma }}[SS{]}_C^{\sigma }{S}^{2\lambda }$(29) ) to demonstrate our method of finding invariants. Suppose we have a function f which satisfies $f(S,[SS{]}_C,\left[\widetilde{SS}{]}_H\right)=0.$ Let $A=\mathrm{\partial }f/\mathrm{\partial }S$, $B=\mathrm{\partial }f/\mathrm{\partial }[SS{]}_C$, and $C=\mathrm{\partial }f/\mathrm{\partial }[\widetilde{SS}{]}_H$. Then, $\frac{\mathrm{d}f}{\mathrm{d}t}=A\frac{\mathrm{d}S}{\mathrm{d}t}+B\frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}+C\frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}=0,$ which by system (24(24) $\begin{array}{rl}\frac{\mathrm{d}S}{\mathrm{d}t}& =-{\beta }_C[SI{]}_C-{\beta }_H[SI{]}_H\\ \frac{\mathrm{d}I}{\mathrm{d}t}& ={\beta }_C[SI{]}_C+{\beta }_H[SI{]}_H-\gamma I\\ \frac{\mathrm{d}[SI{]}_C}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CC}\frac{[SI{]}_C^2}{S}-{\beta }_H{\kappa }_{CH}\frac{\left[SI{]}_C\right[SI{]}_H}{S}-({\beta }_C+\gamma +\delta )[SI{]}_C\\ \frac{\mathrm{d}[SI{]}_H}{\mathrm{d}t}& =-{\beta }_C{\kappa }_{CH}\frac{\left[SI{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{[SI{]}_H^2}{S}-({\beta }_H+\gamma )[SI{]}_H+\eta [\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[\widetilde{SI}{]}_H}{\mathrm{d}t}& ={\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}+{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}\\ & -{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_C}{S}-{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SI}{]}_H\right[SI{]}_H}{S}-(\eta +\gamma )[\widetilde{SI}{]}_H\\ \frac{\mathrm{d}[SS{]}_C}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CC}\frac{\left[SS{]}_C\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{CH}\frac{\left[SS{]}_C\right[SI{]}_H}{S}\\ \frac{\mathrm{d}[\widetilde{SS}{]}_H}{\mathrm{d}t}& =-2{\beta }_C{\kappa }_{CH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_C}{S}-2{\beta }_H{\kappa }_{HH}\frac{\left[\widetilde{SS}{]}_H\right[SI{]}_H}{S}.\end{array}$(24) ) implies $\begin{array}{rl}& {\beta }_C[SI{]}_C\left(A+2B{\kappa }_{CC}\frac{[SS{]}_C}{S}+2C{\kappa }_{CH}\frac{[\widetilde{SS}{]}_H}{S}\right)\\ & +{\beta }_H[SI{]}_H\left(A+2B{\kappa }_{CH}\frac{[SS{]}_C}{S}+2C{\kappa }_{HH}\frac{[\widetilde{SS}{]}_H}{S}\right)=0.\end{array}$ Thus, $df/dt=0$ if the following system of equations is satisfied: $\begin{array}{rl}& A+2B{\kappa }_{CC}\frac{[SS{]}_C}{S}+2C{\kappa }_{CH}\frac{[\widetilde{SS}{]}_H}{S}=0,\\ & A+2B{\kappa }_{CH}\frac{[SS{]}_C}{S}+2C{\kappa }_{HH}\frac{[\widetilde{SS}{]}_H}{S}=0.\end{array}$ For a given C, the solution to this system is: $\begin{array}{rl}A& =-2C\lambda \frac{[\widetilde{SS}{]}_H}{S}\\ B& =-C\sigma \frac{[\widetilde{SS}{]}_H}{[SS{]}_C}.\end{array}$ We choose $C=\mathrm{\partial }f/\mathrm{\partial }[\widetilde{SS}{]}_H=1$ which implies $f=[\widetilde{SS}{]}_H+g([SS{]}_C,S)=0$ for some function g. Therefore, $\frac{\mathrm{\partial }g}{\mathrm{\partial }S}=\frac{\mathrm{\partial }f}{\mathrm{\partial }S}=A=-2\lambda \frac{[\widetilde{SS}{]}_H}{S}=2\lambda \frac{g}{S}$ from which separation of variables gives $g=S^{2\lambda }h\left(\right[SS{]}_C)$ for some function h. Hence, $S^{2\lambda }\frac{\mathrm{d}h}{\mathrm{d}[SS{]}_C}=\frac{\mathrm{\partial }g}{\mathrm{\partial }[SS{]}_C}=\frac{\mathrm{\partial }f}{\mathrm{\partial }[SS{]}_C}=B=\sigma \frac{g}{[SS{]}_C}=\sigma S^{2\lambda }\frac{h}{[SS{]}_C}$ which gives $\frac{\mathrm{d}h}{h}=\sigma \frac{\mathrm{d}[SS{]}_C}{[SS{]}_C}$ from which separation of variables gives $h=K[SS{]}_C^{\sigma }$ for any constant K. Therefore, $\mathrm{d}f/\mathrm{d}t=0$ is satisfied by any f of the form: $f=[\widetilde{SS}{]}_H+g([SS{]}_C,S)=[\widetilde{SS}{]}_H+S^{2\lambda }h(\left[SS{]}_C\right)=[\widetilde{SS}{]}_H+K[SS{]}_C^{\sigma }{S}^{2\lambda }.$ Substituting the initial conditions (25(25) $\begin{array}{rl}S\left(0\right)& ={\alpha }^S,I\left(0\right)={\alpha }^I,R\left(0\right)={\alpha }^R\\ \left[SI{]}_C\right(0)& ={\mu }_C{\alpha }^S{\alpha }^I,\left[SI{]}_H\right(0)=0,[\widetilde{SI}{]}_H\left(0\right)={\mu }_H{\alpha }^S{\alpha }^I\\ \left[SS{]}_C\right(0)& ={\mu }_C({\alpha }^S{)}^2,[SS{]}_H\left(0\right)=0,\left[\widetilde{SS}{]}_H\right(0)={\mu }_H({\alpha }^S{)}^2.\end{array}$(25) ) into f=0 gives $K=-\frac{\left[\widetilde{SS}{]}_H\right(0)}{\left(\right[SS{]}_C\left(0\right){)}^{\sigma }\left(S\right(0){)}^{2\lambda }}=-\frac{{\mu }_H{\alpha }_S^{2(1-\sigma -\lambda )}}{{\mu }_C^{\sigma }}.$ Hence, the invariant is $[\widetilde{SS}{]}_H-\frac{{\mu }_H{\alpha }_S^{2(1-\sigma -\lambda )}}{{\mu }_C^{\sigma }}[SS{]}_C^{\sigma }{S}^{2\lambda }=0,$ i.e. relation (29(29) $[\widetilde{SS}{]}_H=\frac{{\mu }_H({\alpha }^S{)}^{2(1-\sigma -\lambda )}}{{\mu }_C^{\sigma }}[SS{]}_C^{\sigma }{S}^{2\lambda }$(29) ).

The derivations of relations (32(32) $\begin{array}{rl}[SI{]}_C& =\frac{{\beta }_C+\gamma +\delta }{{\beta }_C(1-{\kappa }_{CC})}S-[SS{]}_C+\left(\frac{({\beta }_C+\gamma +\delta )}{{\beta }_C({\kappa }_{CC}-1)}+{\mu }_C({\alpha }^S+{\alpha }^I)\right){\left(\frac{[SS{]}_C}{{\mu }_C}\right)}^{1/2}\end{array}$(32) ) and (33(33) $\begin{array}{rl}[\widetilde{SI}{]}_H& =\frac{({\beta }_H+\gamma )(\eta +\gamma )}{\eta {\beta }_H(1-{\kappa }_{HH})}S-\frac{\eta +\gamma }{\eta }[SI{]}_H-[\widetilde{SS}{]}_H\\ & +\left(\frac{({\beta }_H+\gamma )(\eta +\gamma )}{\eta {\beta }_H({\kappa }_{HH}-1)}+{\mu }_H({\alpha }^S+{\alpha }^I)\right){\left(\frac{[\widetilde{SS}{]}_H}{{\mu }_H}\right)}^{1/2}.\end{array}$(33) ), the additional invariants in the independent layers case, are similar. For (32(32) $\begin{array}{rl}[SI{]}_C& =\frac{{\beta }_C+\gamma +\delta }{{\beta }_C(1-{\kappa }_{CC})}S-[SS{]}_C+\left(\frac{({\beta }_C+\gamma +\delta )}{{\beta }_C({\kappa }_{CC}-1)}+{\mu }_C({\alpha }^S+{\alpha }^I)\right){\left(\frac{[SS{]}_C}{{\mu }_C}\right)}^{1/2}\end{array}$(32) ), we look for an invariant of the form $f(S,[SS{]}_C,\left[SI{]}_C\right)=0$. Likewise, for (33(33) $\begin{array}{rl}[\widetilde{SI}{]}_H& =\frac{({\beta }_H+\gamma )(\eta +\gamma )}{\eta {\beta }_H(1-{\kappa }_{HH})}S-\frac{\eta +\gamma }{\eta }[SI{]}_H-[\widetilde{SS}{]}_H\\ & +\left(\frac{({\beta }_H+\gamma )(\eta +\gamma )}{\eta {\beta }_H({\kappa }_{HH}-1)}+{\mu }_H({\alpha }^S+{\alpha }^I)\right){\left(\frac{[\widetilde{SS}{]}_H}{{\mu }_H}\right)}^{1/2}.\end{array}$(33) ) we consider an invariant of the form $f(S,[\widetilde{SS}{]}_H,[SI{]}_H,[\widetilde{SS}{]}_H)=0$. The analysis follow analogously to above.