{'fraction_non_alphanumeric': 0.02805136164808805, 'fraction_numerical': 0.008903626358120502, 'mean_word_length': 4.741331821127673, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Edge-Based Compartmental Modelling of an SIR Epidemic on a Dual-Layer Static–Dynamic Multiplex Network with Tunable Clustering

The duration, type and structure of connections between individuals in real-world populations play a crucial role in how diseases invade and spread. Here, we incorporate the aforementioned heterogeneities into a model by considering a dual-layer static–dynamic multiplex network. The static network layer affords tunable clustering and describes an individual’s permanent community structure. The dynamic network layer describes the transient connections an individual makes with members of the wider population by imposing constant edge rewiring. We follow the edge-based compartmental modelling approach to derive equations describing the evolution of a susceptible–infected–recovered epidemic spreading through this multiplex network of individuals. We derive the basic reproduction number, measuring the expected number of new infectious cases caused by a single infectious individual in an otherwise susceptible population. We validate model equations by showing convergence to pre-existing edge-based compartmental model equations in limiting cases and by comparison with stochastically simulated epidemics. We explore the effects of altering model parameters and multiplex network attributes on resultant epidemic dynamics. We validate the basic reproduction number by plotting its value against associated final epidemic sizes measured from simulation and predicted by model equations for a number of set-ups. Further, we explore the effect of varying individual model parameters on the basic reproduction number. We conclude with a discussion of the significance and interpretation of the model and its relation to existing research literature. We highlight intrinsic limitations and potential extensions of the present model and outline future research considerations, both experimental and theoretical.


Introduction
The continual design and development of mathematical models describing epidemic processes on large, complex populations improves our understanding of how diseases and individuals behave during an epidemic, and how preventative measures can be implemented for the greater good. With ever-increasing computational power, models can incorporate increasingly complex features, and model predictions may become more valuable. Nonetheless, any model must tread a careful balance between capturing observed real-world complexity and enabling calculations and conclusions to be drawn with ease. The ultimate epidemiological model must therefore incorporate the behavioural and structural features which significantly influence disease dynamics, whilst being analytically tractable.
Social heterogeneity describes the propensity for a social group to be diverse in character or content and is an important determinant when studying the dynamics and control of infectious diseases (Arthur et al. 2017). In a social group, heterogeneity encompasses many descriptive elements, such as variations in individuals' behaviour or in susceptibility across group members. In network theory, social heterogeneity can also describe variations in the types of connections an individual makes. For example, an individual can be connected to other individuals in distinct groups, such as workplace or community groups.
Structured populations with multiple connection types are well described by multiplex networks, where a population of individuals partakes in multiple network layers. Each network layer describes a specific type of interaction between members of the population, and network structure in one layer is allowed overlap with network structure in another layer. A pair of individuals in a multiplex network can share more than one connection. In a multiplex network, an individual is present in every network layer, but may or may not partake in connections in individual network layers.
Existing multiplex modelling studies have shown that single-layer approximations or aggregations of multiplex networks are not accurate enough to describe the epidemic process (Diakonova et al. 2016;Zhuang et al. 2017;Gomez et al. 2013;Cozzo et al. 2013), and further that an epidemic can spread on a multiplex network even if the individual layers are well below their respective epidemic thresholds (Zhao et al. 2014). A global cascades model generalised for multiplex networks was used to show that multiplexes are more vulnerable to global cascades than single-layer networks (Brummitt et al. 2012). These studies highlight the importance of accounting for heterogeneity in connection type by considering multiplex network models.
Another determinant of infectious disease dynamics is heterogeneity in the structural connections between individuals, within a single type of connection. Real-world networks often exhibit community structure, with a high density of connections within communities and a low density of connections between communities. They are also considered to exhibit other structural characteristics such as network transitivity or clustering, described in social network theory as the propensity for an individual to be connected to a friend of a friend (Newman 2003).
Community structure has been shown to affect disease dynamics on single-layered (uniplex) networks, where on average, epidemics occurring on networks with community structure exhibit greater variance in final epidemic size, a greater number of small, local outbreaks that do not develop into epidemics and higher variance in the duration of the epidemic (Salathé and Jones 2010). Network quality functions able to detect community structure in multiplex networks have been developed (Mucha et al. 2010). Further, results such as the large graph limit of a susceptible-infected-recovered (SIR) epidemic process on a dynamic multilayer network, where one network layer represents community links and another represents connections in healthcare settings, have been derived (Jacobsen et al. 2016).
In network models, increased clustering is generally considered to slow an epidemic by increasing the epidemic threshold (Miller 2009). However, this relationship is not always monotonic. Higher clustering in a multiplex study of information propagation led to an increase in the epidemic threshold and a decrease in final epidemic size (Zhuang and Yagan 2016). Increased clustering in a study of Watt's threshold model generalised for a multiplex network comprised of clustered network layers led to a decrease in the probability of a global cascade and its size (Zhuang et al. 2017). However, the authors also discovered a critical threshold for the average degree, above which clustering was shown to facilitate global cascades (Zhuang et al. 2017). A uniplex network study found that simultaneously increasing clustering and the variance of the degree distribution led to an increase in final epidemic size (Volz et al. 2011). Moreover, clustering can lead to correlations where high-degree individuals are more likely to connect with other high-degree individuals. It is clear that the effect of clustering is complex and should be considered in the design of network models.
In epidemiology it is also important to consider heterogeneity across contact duration. In human populations, links between individuals may be long-lasting (persistent), e.g. between an infant child and their caregiver; temporary (transient), e.g. between workplace colleagues; or more short-lived (fleeting), e.g. between strangers coming into close proximity on public transport. In a study using a year's mobile phone data as a proxy for the structure and dynamics of a large social network, researchers found that persistent links tend to be reciprocal and are more common for individuals with low degree and high clustering (Hidalgo and Rodriguez-Sickert 2008). Many network-based studies in the past have considered fully static network structures and hence solely investigate the effects of persistent connections between individuals, see Keeling and Eames (2005) for a review of differing approaches.
Later studies of epidemic processes on networks have incorporated persistent and transient connections into their models by imposing rewiring rules on static networks. Rewiring rules considered include spatially constrained rewiring (Rattana et al. 2014), random link activation and deletion (Taylor et al. 2012;Sélley et al. 2015;Kiss et al. 2012) and temporary link deactivation (Tunc et al. 2013;Shkarayev et al. 2014). On the other hand, epidemic processes with fleeting contact duration can be well described via the mass action model, which assumes all pairs of individuals contact one another at the same rate, the mean-field social heterogeneity model (also known as the degree-based mean-field model), which generalises the mass action model by allowing for variations in contact rate across the population, and the dynamic fixedand dynamic variable-degree models, where edges are swapped at a given rate, or edges are broken and created at given rates, respectively (Miller and Volz 2013;Miller et al. 2012).
Here, we suppose that static and dynamic connections coexist in any complex population. We aim to derive a network model describing an SIR epidemic process spreading through a population where each individual has two types of connections: persistent links to individuals in their household, constituting a static network layer with community structure, and transient connections to strangers in the wider population, where all such edges rewire at a constant rate, constituting a dynamic network layer with conserved degrees.
In what follows, we utilise the edge-based compartmental modelling (EBCM) approach (Volz 2008;Miller 2011Miller , 2014Miller et al. 2012), deriving equations which describe the time evolution of classical quantities of interest, where the underlying dual-layered static-dynamic network has heterogeneity in contact type, contact duration, and contact structure. We derive the associated basic reproduction number R 0 , following the next-generation matrix approach (Diekmann et al. 2009). We describe the implementation of the EBCM model and of statistically correct Gillespie simulations of the epidemic process (Gillespie 1976). The new model is validated, firstly by showing that collapsing either the static or dynamic network layers leads model equations to converge to existing equivalent model equations, and secondly by comparing the dynamics predicted by model equations to those from exact simulations. We explore how various combinations of model parameters and network layers influence global dynamics, uncover behavioural regimes that the model can achieve for specific combinations of infection and rewiring rates, and show that our derived R 0 behaves as expected. The paper concludes with a discussion of potential implications of the work as well as possible extensions.

Methods
Our solutions are based on the class of undirected random graphs (networks). Each node is a member of a random number of static lines (2-vertex cliques), static triangles (3-vertex cliques) and dynamic lines (2-vertex cliques). The probability that a node has s static line stubs, t static triangle corners and d dynamic line stubs is described by the probability mass function p s,t,d . The model captures network structure using the probability generating function (PGF) When differentiating PGF (1), we use superscripts such that g (x) denotes the first (partial) derivative of g with respect to x and g (y,y) denotes the second (partial) derivative of g with respect to y. Equation (1) can be used to calculate useful properties of the multiplex network. For example, M s , the expected number of static line stubs that belong to a randomly selected individual, M t , the expected number of static triangle corners that belong to a randomly selected individual, and M d , the expected number of dynamic line stubs that belong to a randomly selected individual, are calculated as follows: We consider a basic SIR compartmental model. Infections occur across edges on the static network layer at a constant rate β s , whilst infections occur across edges on the dynamic network layer at a constant rate β d . Infected individuals recover at a constant rate γ . Once recovered, a node cannot be reinfected, and can no longer transmit infection to its neighbours. A comprehensive list of model variables and parameters is given in Table 1.

Edge-Based Compartmental Model Derivation
We follow the edge-based compartmental modelling approach by considering the fate of a randomly selected test node u, which is prevented from transmitting infection. This assumption is a useful tool that eliminates conditional probability arguments that would need to be considered otherwise (Miller et al. 2012). It does not introduce any approximation. At time zero, infection is introduced to a fraction ρ of the population chosen uniformly at random, comprising the initial condition of the system. We assume that the test node u is a member of s static line stubs, t static triangle corners and d dynamic line stubs. Then the probability that u is susceptible is (1 − ρ)θ s 2 θ t 3 θ d 4 , where θ 2 is the probability that a random line (2-clique) on the static network layer has not transmitted infection to the test node, θ 3 is the probability that neither of the other nodes in a random triangle on the static network layer has transmitted infection to the test node, and θ 4 is the probability that a random stub connected to u on the dynamic network layer has never been involved in transmitting infection to the test node. Assuming we are able to calculate θ 2 , θ 3 and θ 4 as functions of time, we are able to calculate the proportion of susceptible individuals S as a function of time. Given S(t), we use I (t) = 1 − S(t) − R(t) andṘ(t) = γ I (t) to calculate I (t) and R(t), completing the system.

Considering Â 2
We divide θ 2 into φ S , φ I and φ R , the probabilities that a random neighbour along a line on the static network layer has not transmitted infection to u, and is susceptible, infected or recovered, respectively. The probability the neighbour has not transmitted infection to u is θ 2 = φ S +φ I +φ R , and (1−θ 2 ) is the probability that it has transmitted infection to u. The fluxes between these quantities are shown in Fig. 1. The fluxes from φ I to φ R and from φ I to (1 − θ 2 ) are proportional to one another. Both φ R and (1 − θ 2 ) are equal to zero at time zero since we assume that no infection or recovery events A survivor function for remaining susceptible for some time t ≥ 0, given that the individual in question is a member of a single static line θ 3 (t) A survivor function for remaining susceptible for some time t ≥ 0, given that the individual in question is a member of a single triangle corner θ 4 (t) A survivor function for remaining susceptible for some time t ≥ 0, given that the individual in question is a member of a single dynamic line The probabilities that a neighbour of u along a static line is susceptible, infectious or recovered, and has not transmitted infection to u by time t ≥ 0 φ XY (t) The probability that two neighbours of u in a triangle are in states X and Y ∈ {S, I , R} and have not transmitted infection to u by time t ≥ 0 The rate (at time t ≥ 0) at which a random triangle neighbour v of u is infected from outside the triangle, given that v was susceptible The rate (at time t ≥ 0) at which a random dynamic line stub neighbour v of u becomes infected from outside the dynamic line joining u and v, given that v was susceptible The probabilities (at time t ≥ 0) that a random dynamic stub belonging to u has never been involved in transmitting infection to u and is currently connected to a susceptible, infected or recovered individual, respectively π S (t), π I (t), π R (t) The probabilities (at time t ≥ 0) that a randomly chosen dynamic stub belongs to a susceptible, infected or recovered individual, respectively Many definitions refer to the test node u, which is selected at random from the population and modified so that it cannot transmit infection, but can itself become infected Fig. 1 Flow diagram for the flux of a static line partner through different states. The flux between the probabilities that the test node u is connected by a line (2-clique) on the static network layer to a node v that has not transmitted infection to u and is susceptible (φ S ), infectious (φ I ) or recovered (φ R ), and the probability that v has transmitted infection to u, equal to (1 − θ 2 ) can occur prior to time zero. By integrating the relation dφ R dt = γ β s d(1−θ 2 ) dt , and using the initial condition φ R (0) = (1 − θ 2 (0)) = 0, we find the relation Next, we must calculate an expression for φ S . Consider the number of static line stubs attached to an individual that we reach by following a randomly chosen static line. Similarly, consider the number of static triangle corners attached to an individual reached by following a randomly chosen static triangle edge, and the number of dynamic line stubs attached to an individual we reach by following a randomly chosen dynamic line. Following edges in this way means we are more likely to arrive at an individual with a higher degree, in direct proportion to that individual's degree (Meyers et al. 2006). The random number of such lines and triangle corners is described by the excess degree distribution, and we calculate the associated probability density functions for each edge type as follows. Denote q s−1,t,d ∝ sp s,t,d as the probability of there being (s − 1) static line stubs, t triangle corners and d dynamic line stubs connected to a susceptible node that we reach by following a static line, not counting the line by which we arrived. Similarly, denote r s,t−1,d ∝ t p s,t,d as the probability that if we follow a triangle edge to a susceptible node, there are s static line stubs, (t − 1) triangle corners and d dynamic line stubs connected to that node, not counting the triangle edge by which we arrived, and w s,t,d−1 ∝ dp s,t,d as the probability that if we follow a dynamic edge to a susceptible node, there are s static line stubs, t triangle corners and (d − 1) dynamic line stubs connected to that node, not counting the dynamic edge by which we arrived.
From above, we note that the probability that there are s static line stubs, t triangle corners and d dynamic line stubs attached to a random neighbour of u across a static line (not counting the line it was reached across) is q s−1,t,d ∝ sp s,t,d . A neighbour reached by following a static line connected to u is susceptible with probability (1−ρ)θ s−1 2 θ t 3 θ d by M s = g (x) (1, 1, 1), the expected number of static lines a randomly selected node belongs to. We find From the original definition of θ 2 we have We are now able to calculate an expression for θ 2 using Eqs.

Considering Â 3
Since θ 3 denotes the probability that neither of the other nodes in a triangle has transmitted infection to the test node, we must divide θ 3 into six quantities φ SS , φ SI , φ SR , φ II , φ IR and φ RR in order to consider all possible disease status combinations for two individuals. For example, φ SI denotes the probability that one triangle neighbour of u is susceptible, whilst the other is infectious, and neither has transmitted infection to u. The flux between the various compartments can be seen in Fig. 2. There is no simple relation between φ RR and θ 3 , so we take a different approach than before. We start withθ 3 , which satisfieṡ To calculate elements in the right-hand side of (6), we must first obtain an expression for φ SS , the probability that both neighbours in a triangle are still susceptible. Under the assumption that no transmission events have occurred in the triangle, the probability that a single triangle neighbour of u is susceptible is where M t is the expected number of static triangle corners belonging to a randomly chosen individual. Since we require both triangle neighbours of u to be susceptible, we have We choose A to denote the rate at which a single triangle neighbour of u becomes infected from outside the triangle. From Fig. 2 we know that dφ SS dt = −2 Aφ SS , which

Fig. 2
Flow diagram for the flux of two triangle neighbours through different states. The flux between the probabilities that the test node u is connected in a triangle to two nodes in all possible disease status configurations, where neither triangle neighbour has transmitted infection to u, as well as the probability (1 − θ 3 ) that a node v = u in the triangle has transmitted infection to the test node u implies A = − dφ SS dt /2φ SS . To arrive at an explicit formula for A, we begin by calculating dφ SS dt via the chain rule: We know that g (y) (1, 1, 1) = 0, since g (y) (1, 1, 1) = s,t,d t p s,t,d ∈ R. Hence (1, 1, 1) .
Using A = − dφ SS dt /2φ SS and some simplification, we find an explicit formula for A: Now we are ready to calculate equations for φ SI , φ II and φ IR . We also require φ SR , but do not require a formula for φ RR . Using the flow diagram in Fig. 2, we havė

Considering Â 4
To take into account the dynamic rewiring of edges, we introduce θ 4 = ψ S + ψ I + ψ R , where ψ I denotes the probability that a random dynamic stub belonging to the test node u has never been involved in transmitting infection to u, and is currently connected to an infectious node. Other important assumptions with respect to dynamic-edge rewiring are the following: we assume that when one partnership ends, a new partnership forms immediately, neglecting any between-partner period, and we assume that edges break at rate η. The flux between the various compartments of interest can be seen in Fig. 3. Previously, φ S (which corresponds to ψ S in this subsection) was calculated explicitly as the probability that the neighbour is susceptible. With dynamic-edge rewiring, an edge that previously transmitted infection may later become connected to a susceptible node, so the previous calculation of φ S does not apply here. To find ψ S , we need to calculate the probability that a newly formed edge connects to a susceptible, infectious or recovered individual. We call these probabilities π S , π I and π R and note that they are equivalent to the probabilities that a randomly chosen dynamic stub belongs to a node in each disease compartment. The flux between these probabilities can be seen in Fig. 4.
First, we calculate the values π S , π I and π R , beginning with π S . If we select a dynamic stub at random, the probability that it belongs to an individual partaking in s static lines, t triangles and d dynamic stubs is dp is the expected number of dynamic edges that a random individual belongs to. At time zero, infection is introduced at random to a proportion ρ of the population. Thus the probability of any node being susceptible at time zero is (1 − ρ). The probability of a node with degree (s, t, d) being susceptible after some time, given that it was susceptible at time zero, is θ s Flow diagram for the flux of a dynamic-edge partner through different states. The flux between the probabilities θ 4 = ψ S + ψ I + ψ R that a random stub currently connected to u on the dynamic network layer has never been involved in transmitting infection to u. Note that the compartment denoted ηθ 4 is not a compartment in the typical sense. When edges break (at rate η) in the model, moving into 'compartment' ηθ 4 , new edges are formed immediately without delay, moving straight back into compartments ψ S , ψ I or ψ R . π S , π I and π R denote the probabilities that a randomly chosen dynamic stub belongs to a susceptible, infected or recovered node, respectively π S π I π R γπ I Fig. 4 Flow diagram for the flux of a dynamic line stub through different states. The flux between π S , π I and π R , the probabilities that a randomly chosen dynamic stub belongs to a susceptible, infected or recovered node, respectively with the summation taken over all degree possibilities described by the probability mass function p s,t,d . Stubs belonging to infected nodes become stubs belonging to recovered nodes at rate γ ; hence,π R = γ π I , and π I = 1 − π S − π R . The equation for π S can be condensed using PGF (1), so we have To complete the system we need to calculate the flux Bψ S from ψ S to ψ I by solving a differential equation for ψ S . B describes the rate at which a susceptible dynamicedge neighbour v of u becomes infected from outside the dynamic edge joining u and v. Consider a random test node u and a random dynamic-edge neighbour v of u, at some time t. Let ζ denote the probability that the two stubs joining u and v have not previously been involved in transmitting infection to u or to v, prior to the u − v edge forming. The probability that v is susceptible and that u's stub has not previously transmitted to u is , where s is the number of static lines v partakes in, t is the number of triangles v partakes in, and d is the dynamic line stub degree of v. Since we do not know the values (s, t, d) for v, we must consider all possible combinations of degrees. The probability of a randomly chosen dynamic stub belonging to a node with degree (s, t, d) is dp s,t,d /g (z) (1, 1, 1). We conclude that To calculate the derivative of ψ S , we first consider the derivative of ζ . This is given by subtracting the rate at which such edges break, ηζ , from the rate at which such edges form, ηθ 2 4 (one θ 4 for u's stub and one for v's stub). We havė We have an expression forζ , so the derivative of ψ S can be found via the chain rule: Fig. 3 we haveψ S = ηθ 4 π S − ηψ S − Bψ S , so we calculate the flux between compartments ψ S and ψ I using the rate The ψ S to ψ I flux is the product of ψ S , the probability that a random dynamic stub has not transmitted infection to the test node u and is currently connected to a susceptible node, with rate B, the rate that a neighbouring susceptible node v becomes infected from outside the dynamic edge, given that the stub has not transmitted and connects u to a susceptible node. Following the flow diagram in Fig. 3, we have the differential equationsθ

Population-Level Equations
We began the EBCM derivation by considering the probability of a randomly selected test node u (which is prevented from transmitting infection) being susceptible as θ s 2 θ t 3 θ d 4 , given that the node has degree (s, t, d). Since we have calculated formulae for θ 2 , θ 3 and θ 4 , we can derive population-level equations describing the proportion of the population in each disease compartment at each point in time: Equations (1)-(23) form a complete system describing an SIR epidemic spreading across a dual-layer multiplex network consisting of a static network layer constructed from line stubs and triangle corners and a dynamic network layer constructed from line stubs only, where edges rewire and degrees are conserved.

Deriving the Basic Reproduction Number R 0
The basic reproduction number R 0 is defined as the average number of infections caused by a single infectious individual, early in an epidemic process, in an otherwise susceptible population. In the model, a multiplex network structure is generated using three distinct edge distributions (static line stubs, static triangle corners and dynamic line stubs). To compute R 0 we must consider the average number of infections caused across each type of edge, whilst also considering the type of edge that the infection was originally received across. With 3 edge types, this constitutes 9 values, grouped together to form the next-generation matrix where matrix element G i j describes the average number of infections caused across edges of type j, where the infector received infection across an edge of type i. Following the next-generation matrix approach (Diekmann et al. 2009), the value of R 0 is found via the leading eigenvalue of the matrix G, or equivalently, the eigenvalue with greatest magnitude. We note that the matrix G defined here is the transpose of the next-generation matrix as defined in Diekmann et al. (2009). However, this discrepancy does not affect the eigenvalues or therefore the value R 0 . To find R 0 , we begin by deriving expressions for values in the first column of G. Firstly, consider the non-diagonal matrix entries G ts and G ds . We want to compute the expected number of infection events occurring across static lines, when individuals contracted infection across a triangle edge or a dynamic line. In both cases, we require the expected static line stub degree, multiplied by the expected number of infections caused across a single static line attached to the infectious individual. Say the expected static line stub degree is denoted k s . Now we require the expected number of infections caused across a single static edge attached to an infectious individual, in an otherwise susceptible population. A single static edge joining a susceptible and an infectious individual, in an otherwise susceptible population, has two event possibilities: a single recovery, or a single infection. Denote X as the random variable describing the number of infection events occurring across a single static line joining a susceptible to an infectious individual, in an otherwise susceptible population. Using the expectation formula, and since there can only be zero or one infection events occurring across such an edge, we find the expected number of infections across a static line joining a susceptible to an infectious individual simply as P(X = 1). The probability of a single infection occurring across such a static edge, prior to any recovery, is β s β s +γ . Thus we can say that G ts = G ds = k s β s β s +γ . Finally, we calculate an expression for the diagonal matrix element G ss , by multiplying the expected excess static line stub degree, denoted s , by the expected number of infections caused across a single static line joining a susceptible individual to an infectious individual in an otherwise susceptible population. Following the same argument for G ts and G ds , we compute the expected number of infection events for G ss as β s β s +γ , and we obtain G ss = s β s β s +γ . Next we derive expressions for the values G st , G tt and G dt in the second column of the matrix G. We firstly consider the non-diagonal elements G st and G dt . Both G st and G dt are calculated by multiplying the expected triangle corner degree, denoted k t , by the expected number of infection events caused within a single triangle attached to an infectious node in an otherwise susceptible population. In a single triangle comprised of two susceptible individuals attached to an infectious individual, there are a finite number of infection event possibilities: either no further infections occur (the infectious individual recovers), one infection event occurs, or two infection events occur. Define Y as the random variable describing the number of infection events within such a triangle. Using the expectation formula, we find the expected number of infection events within a triangle comprised of two susceptible individuals and an infective, in an otherwise susceptible population, as P(Y = 1)+2·P(Y = 2). To continue, we must compute the probabilities P(Y = 1) and P(Y = 2) explicitly. P(Y = 1) describes the probability that the original infective infects one out of two triangle neighbours. In this case, either one of the two susceptible neighbours can become infectious, and both infectious triangle members must then recover, so that it is impossible for any more than one infection event to occur. In a triangle comprised of a single infective and two susceptible nodes, there are four distinct nodal orders in which a single infection event is followed by the recovery of both infectious nodes. We find Considering P(Y = 2) is more complex, as there are two distinct ways in which two infection events can occur in a triangle between an infective and two susceptible individuals. Firstly, the original infective can infect both of its triangle neighbours consecutively, prior to any recovery events. The probability of both triangle infection events occurring in succession is given by Secondly, the original infective can cause two triangle infections via three consecutive events. In this case, the originally infectious triangle member firstly infects one susceptible triangle neighbour at rate 2β s 2β s +γ . The triangle is now comprised of two infectious individuals attached to a single susceptible individual. The second event to occur is a recovery of either the original infector or its first infectee, occurring at rate 2γ 2β s +2γ = γ β s +γ . The triangle is now comprised of a susceptible, an infective, and a recovered individual, in an otherwise susceptible population. Following the recovery event, the final event is an infection of the remaining susceptible triangle member, occurring at rate β s β s +γ . The probability of all three events occurring in succession is thus 2β s 2β s +γ γ β s +γ β s β s +γ . In the latter case of an infection, followed by a recovery, followed by another infection within a triangle originally composed of an infective and two susceptible individuals in an otherwise susceptible population, the original infector may not be directly involved in every single infection event. However, for the purposes of deriving R 0 , we say that the original infector caused these infections, regardless of the order in which triangle members recover and infect one another.
Since there are two distinct ways in which two infections can take place within a triangle comprised of an infective and two susceptible individuals, we take the sum of both individual probabilities to obtain P(Y = 2): We find the expected number of infection events within a triangle comprised of two susceptible individuals and an infective, in an otherwise susceptible population, as We conclude by deriving elements from the third column of G, starting with nondiagonal matrix elements G sd and G td . In both cases, we multiply the expected dynamic line stub degree, denoted k d , by the expected number of infection events occurring across a single dynamic line stub attached to an infectious individual, in an otherwise susceptible population.
The probability of a dynamic stub attached to an infective in an otherwise susceptible population transmitting infection at least once is β d β d +γ . If such an infection occurs, the I-S pairing becomes an I-I pairing with a dynamic edge joining the two individuals. The probability of a dynamic I-I edge rewiring, prior to any recovery event, is η η+γ . We can assume that any I-I edge rewires to become an I-S edge in the limit of large population size, since we are early on in an epidemic process, and we began with an otherwise susceptible population. The probability that an infectious dynamic stub infects its new susceptible neighbour is β d β d +γ . This rewiring and infecting process can occur an arbitrary number of times in the model. The expected number of infections of this type can be calculated by taking the sum by the geometric series, and where r is defined as , the probability of an infectious individual's dynamic-edge rewiring, followed immediately by its dynamic stub infecting the new (susceptible) neighbour across the rewired edge. We obtain the matrix values G sd = k d β d β d +γ 1 1−r = G td . Finally, we compute G dd , defined as the expected number of infections caused across dynamic edges, where the infector received infection across a dynamic edge itself. Firstly, consider the single dynamic I-I edge which originally infected our individual. The probability of the edge rewiring, leaving our infective in an I-S dynamic-edge pairing, is η η+γ . The probability of the infectious dynamic stub infecting the new susceptible neighbour is β d β d +γ . Thus the probability that the dynamic stub which originally contracted infection infects ≥ n individuals is r n , where . We compute the expected number of infections of this type by taking the sum of r n for n = 1:∞ ∞ n=1 r n = r 1 − r , by the geometric series. Now consider the remaining dynamic edges associated with our infectious individual. We require the expected number of infections caused by a single edge of this type. Using the same argument as for G sd and G td , we find the expected number of infections caused by one dynamic edge attached to our infectious individual as β d where d is the expected excess dynamic line stub degree.
In detail, the next-generation matrix G takes the form ⎛ where k s , k t and k d denote the expected static line stub, static triangle corner and dynamic line stub degrees, s , t and d denote the expected excess static line stub, static triangle corner and dynamic line stub degrees, and r = ηβ d (η+γ )(β d +γ ) . The basic reproduction number R 0 is the eigenvalue of next-generation matrix (24) with greatest magnitude.

Model Implementation
A variable-order stiff differential equation solver (ode15s in the MATLAB environment) was used to solve all relevant systems of equations. Initial conditions were specified, consisting of appropriate degree distributions and parameters for each edgebased compartmental model type, and of a user-specified end time for the computation.
Solutions to Eqs. (1)-(23) were found using both interdependent and independent distributions for the three edge types. For interdependent distributions, a single probability distribution governed the distribution of pairs of edge stubs, and additional model parameters ( p s + p t + p d ) ≡ 1 were used to distribute each pair of stubs into: two static line stubs (with probability p s ), a single static triangle corner (with probability p t ) or two dynamic line stubs (with probability p d ). In such cases we used a negative binomial distribution for pairs of edge stubs with parameters p and r describing the probability of success in a single trial and the number of trial successes, respectively, where the distribution itself is generated by g nb (x; r , p) = ( p 1−(1− p)x ) r and models the number of failures before a specified number of successes is reached in a series of identical, independent Bernoulli trials. We also utilised a discrete homogeneous distribution for pairs of edge stubs where all individuals had identical degree.
For independent distributions, we used three separate binomial distributions for the number of static line stubs, static triangle corners and dynamic line stubs.

Simulation Implementation
To test the validity of solutions to Eqs. (1)-(23), found in the MATLAB environment, Gillespie simulations (Gillespie 1976) were implemented to produce statistically correct trajectories of SIR epidemic processes occurring on equivalent static-dynamic multiplex networks. Prior to each simulation, static and 'dynamic' adjacency matrices were generated according to a configuration model approach, described as follows: for a population of N individuals, three vectors of length N were generated to record the number of static line stubs, static triangle corners and dynamic line stubs associated with each individual, according to user-specified degree distributions provided to the script. The script ensured that the total number of static line stubs was even, the total number of dynamic line stubs was even, and that the total number of static triangle corners was a multiple of three.
Firstly, the static network layer was generated using vectors containing the number of static line stubs and triangle corners each individual partook in. Pairs of static line stubs and triples of static triangle corners were selected at random. Provided potential static lines and triangles did not generate self-loops (where an individual is joined to itself with an edge) or double edges (where an edge exists more than once within the static network layer), they were added to the static adjacency matrix. The unmatched static line stubs and static triangle corners lists were updated, and the process continued until all static line stubs and triangle corners were successfully matched.
Secondly, the initial structure of the dynamic network layer was generated using the vector storing the number of dynamic line stubs each individual partook in. Pairs of dynamic line stubs were selected at random. Provided a potential dynamic edge did not generate a self-loop or a double edge within the dynamic network layer, it was added to the dynamic adjacency matrix. Successfully paired dynamic stubs were removed from the unmatched stubs list, and the process continued until all dynamic line stubs were successfully matched.
The nature of this configuration model approach meant that the wiring processes for the static and dynamic network layers may have had to be restarted multiple times in order to achieve final network structures. Once all static line stubs, static triangle corners and dynamic line stubs had been wired up, the configuration process was complete. Although the script prevented double edges from occurring within each network layer, it was possible for double edges to occur across the network layers, i.e. for two individuals to share both a static and a dynamic connection simultaneously.
Given static and dynamic adjacency matrices describing the multiplex network structure, simulated epidemic processes were implemented. In each Gillespie simulation, ρ N initially infectious individuals were selected at random from the population. At each time step, a vector of length (N + 1) described the state transition rate (infection or recovery) for all N individuals, followed by a single edge swapping rate, ηM 2 , where M := total number of edges in the dynamic network layer. Inter-event times followed an exponential distribution with scale parameter 1 R , where R := the sum of the rates vector at the current time step. Each event occurring was either an infection, a recovery or an edge swap. Uniformly distributed random numbers were generated at each time step to determine the next event to occur. When an edge swap event occurred, the script selected two dynamic edges at random, ensuring that all four nodes involved in these edges were unique. The script also ensured that the proposed new dynamic edges did not already exist within the dynamic network layer. Given these conditions, an edge swap occurred and the Gillespie process continued. The process terminated once the user-specified end time was reached.

Results
In what follows, we assess the validity of Eqs.
(1)-(23) and of the basic reproduction number R 0 , obtained via next-generation matrix (24). We firstly consider two extreme cases of the multiplex model: when either the static or the dynamic network layers are negligible (close to zero). In such cases, we show that predictions made by Eqs.
(1)-(23) resolve to predictions made by existing uniplex EBCM equations. When the full multiplex model is considered, with static and dynamic network elements present, there exists no basis for comparison other than generating exact simulations of the epidemic process. To this end, we utilise Gillespie simulations to demonstrate the validity of Eqs.
(1)-(23) in predicting the epidemic process for a number of multiplex network configurations. By solely considering the predictions of Eqs.
(1)-(23), we explore the consequences of varying individual model parameters and of considering various combinations of model parameters ( p s + p t + p d ) ≡ 1, governing the contributions of each edge type. Further, we explore the contributions of each edge type and how the resulting final epidemic size is altered within a systematic consideration of combinations of model parameters β s , β d and η. Finally, we test the performance of the derived basic reproduction number R 0 in predicting the outcome of an epidemic and we explore variations in the value of R 0 and the associated final epidemic size predicted by Eqs.
(1)-(23) when altering the rate of rewiring, the extent of clustering and the average degree in the multiplex model.

Model Without Dynamic Layer
When the dynamic component of the dual-layer static-dynamic multiplex is removed, the model reduces to describe an SIR epidemic on a static uniplex network generated by lines and triangles. Biologically speaking, this reduced model tracks the epidemic as it spreads across persistent connections in a population with community structure. The EBCM approach has been followed to derive equations describing an SIR epidemic on such a network (Volz et al. 2011).  close to zero, we were able to test the multiplex model's convergence (Fig. 5). Excellent agreement was observed between multiplex model equations where dynamic network elements are negligible, uniplex model equations (Volz et al. 2011) and Gillespie simulated epidemics on equivalent multiplex networks, for a number of scenarios with varying forces of infection.

Model Without Static Layer
When the static component of the dual-layer static-dynamic multiplex is removed, the model describes an SIR epidemic on a dynamic uniplex network generated by lines, where edges rewire at constant rate η and degrees are conserved. Biologically, the reduced model describes an epidemic spreading through a population where connections between pairs of individuals are temporary, but the number of connections an individual partakes in remains fixed. The EBCM approach was followed to derive equations describing an SIR epidemic spreading on such a network in Miller et al. (2012).
Excellent agreement was observed between predictions made by Eqs.
(1)-(23) when static network elements are close to zero, output from the dynamic fixed-degree derivation in Miller et al. (2012) and Gillespie simulations describing the SIR epidemic and edge rewire processes occurring simultaneously on equivalent multiplex networks, for a number of set-ups with varying forces of infection (Fig. 6).

Model Validation by Comparison with Simulation
We have observed excellent agreement between multiplex model predictions, uniplex model predictions and Gillespie simulated epidemics in extreme cases where either static or dynamic network elements are negligible (Figs. 5, 6). When multiplex network elements are non-negligible, static and dynamic network layers coexist in the model. In such cases, Gillespie simulated epidemics become the sole basis for assessing the validity of multiplex model equations (1)-(23).
A number of comparisons have been made between multiplex model predictions and Gillespie simulations when static and dynamic network elements coexist (Figs. 7,8). Excellent agreement was observed for a number of comparisons with various average degrees (imposed via negative binomial parameters p and r , describing the distribution governing pairs of edge stubs) and various levels of clustering (imposed by varying parameter p t with the constraint ( p s + p t + p d ) ≡ 1) (Fig. 7). Excellent agreement was also observed for a number of comparisons with various combinations of the multiplex model's infection parameters β s and β d (Fig. 8).

A Brief Exploration of Parameter Spaces
Having observed excellent agreement between simulated epidemic processes and equivalent predictions made by multiplex model equations, we investigated the effects of varying single parameters on the dynamics of epidemics predicted by Eqs. (1)-(23). In total, 9 individual model parameters were varied systematically, whilst all (or the majority of) other parameters were held constant (Fig. 9). Across all parameters being varied, an identical baseline parameter set was utilised, with the resulting prediction made by Eqs. (1)-(23) plotted in black to enable ease of comparison between different parameter scenarios.
This brief exploration highlights the effect that increasing or decreasing a single parameter has on the global dynamics of an SIR epidemic spreading across a dual- layer static-dynamic multiplex. Larger values of p, where p describes the probability of success in a single Bernoulli trial, generate a negative binomial distribution with smaller average degree and a reduction in variance, slowing the epidemic's spread. Larger values of r , where r denotes the number of successful Bernoulli trials that must be reached before the experiment is stopped, led the epidemic to spread more rapidly due to an increase in average degree and variance of the negative binomial distribution for pairs of edge stubs. Varying the rewiring rate η led to less pronounced differences, where larger values of η led to a slight increase in the speed at which the epidemic spread through the population. Increasing a single infection parameter β s or β d leads to an increase in the rate of epidemic spread. Altering the parameter ρ means changing the number of individuals who are infectious at the start of an epidemic process. Increasing the value of ρ leads to changes in the shape of the curve I (t), describing the prevalence of infection at time t, and to the epidemic process finishing sooner. Altering the values p s , p t and p d , with the constraint ( p s + p t + p d ) ≡ 1, demonstrates the range of dynamics that can be achieved using a fixed distribution for pairs of edge stubs with additional parameters to distribute edge pairs into three edge types. Baseline infection parameters are used across all three panels, thus β s = 0.05 < 0.2 = β d , meaning that an increase in the proportion of dynamic edges leads to an increase in the speed of the   In the bottom row of panels, parameters p s , p t and p d are being varied. Since the model has the constraint ( p s + p t + p d ) ≡ 1, we alter the triplet values in each panel in the following way. Assume we are varying the parameter p s . If the new p s is larger than the baseline p s , we subtract 1 2 the difference from the remaining baseline parameters p t and p d . Conversely, if the new p s is smaller than the baseline p s , 1 epidemic, whilst any increase in the proportion of static edges leads to a decrease in the rate of epidemic spread.

Contribution of Network Layers via (p s + p t + p d ) ≡ 1
When degree distributions are interdependent, the parameters ( p s + p t + p d ) ≡ 1 afford the ability to investigate the effects on epidemic dynamics of altering the proportion of edges of each type. Previously, we observed changes in the dynamics of I (t), caused by altering the contributions of each edge type (Fig. 9), where β d > β s , rewiring was slow, and pairs of edge stubs were governed by a negative binomial distribution.
In this multiplex setting, increasing the force of infection on one network layer effectively reduces the force of infection on remaining network layers. Thus the value of parameters β s , β d and η, and the ratios between them, bias the effect of varying model parameters p s , p t and p d . To take this into account, we allowed parameters β s , β d and η to take three distinct values (specifically β s ∈ [0.55, 0.6, 0.65], β d ∈ [ β s 2 , β s , 2β s ] and η ∈ [0.01, 1, 100]), and we considered all 27 combinations of their values, before varying the contributions of each edge type and recording the final epidemic size predicted by Eqs. (1)-(23) in each case (Fig. 10). This approach enabled isolation of the effects of changing single infection or rewiring parameters and exploration of the contributions made by various combinations of edge proportions p s , p t and p d in distinct parameter settings.
Increasing the proportion of triangle corners via p t consistently led to decreases in final epidemic size, suggesting that clustering slows the epidemic process regardless of the choice of parameters β s , β d and η (Fig. 10). Generally, increasing the value of η resulted in an increase in final epidemic size when comparing identical edge contributions. Likewise, increasing the value of infection parameters β s or β d led to an increase in final epidemic size. Dependant on the combination of parameters β s , β d and η, different behavioural regimes emerge, indicated by the orientation of colours and the direction in which they change in individual panels. We observe that a single edge proportion can have a more or less dominant effect on the outcome, dependent on the particular parameter set. For example, when η = 0.01 and β s = 0.55 = β d , changing the proportion of dynamic edges p d has little effect on the final epidemic size. However, when η = 100, β s = 0.65 and β d = 1.3, altering the parameter p d leads to more extreme changes in final epidemic size, a result of β d dominating β s and an increased rate of dynamic-edge rewiring.

Validation of Basic Reproduction Number R 0
Next-generation matrix G (24) and the value R 0 can be validated by testing to see if the final epidemic size is disturbed as R 0 exceeds the epidemic threshold (R 0 = 1). When the basic reproduction number is subthreshold (R 0 < 1), the associated epidemic process is expected to 'die-out'. However, when R 0 > 1 the epidemic is expected to take hold and spread within a population.
For a number of set-ups, we recorded the final epidemic size predicted by Eqs.
(1)-(23), the final epidemic size of a single Gillespie simulation of the same process and the associated R 0 value (Fig. 11). To obtain a suitable range of R 0 values we systematically increased β s = β d from subthreshold values, whilst all other parameters were held constant. Independent binomial distributions were used for static line stubs, static triangle corners and dynamic line stubs. In Gillespie simulations where R 0 > 1, we imposed an additional constraint requiring the number of infectives to reach at least ten times the initial number of infected individuals; otherwise, a new Gillespie simulation was implemented. As R 0 exceeded the epidemic threshold, the final epidemic size predicted by model equations (1)-(23) and from individual simulations increased rapidly, suggesting the derivation of the next-generation matrix G and associated R 0 is rigorous. We plotted R 0 and the associated final epidemic size predicted by Eqs.
(1)-(23) for a number of scenarios to investigate the impact on their values of varying specific multiplex network attributes (rewiring, clustering and average degree) and to explore the relationship between R 0 and final epidemic size (Fig. 12). Varying the rewiring rate η demonstrates that R 0 and the associated final epidemic size increase with the value of η. Varying η can also move the system below or above the epidemic threshold R 0 = 1. However, there is a limit to this relationship; as η increases above 20, the changes in R 0 and final epidemic size are negligible. We have seen previously that larger values of p t result in smaller final epidemic sizes, suggesting that increased clustering slows epidemic processes on multiplex networks (Fig. 10). Here, we find that increasing p t leads to decreases in both R 0 and the associated final epidemic size (Fig. 12). The relationship between p t and final epidemic size appears to be linear. For smaller p t the curve with R 0 appears to be linear, but as p t tends towards its maximal value, the reduction in R 0 increases. An increase in average degree k , where pairs of edge stubs follow a negative binomial distribution, led to increases in R 0 and final epidemic size (Fig. 12). The relationship between k (negative binomial) and R 0 appears to be linear. However, the relationship between k and final epidemic size differs. The final epidemic size increases at a faster rate above some critical average degree, say k = 12. A similar pattern emerges in the relationship between the average degree, R 0 and final epidemic size when pairs of edge stubs follow a discrete homogeneous distribution. This is not surprising, as we saw previously that the relationship between R 0 and final epidemic size is nonlinear (Fig. 11). However, these results show that small average degrees make it hard for the epidemic to take hold in the population. Potentially, this is a result of the multiplex network becoming divided into more than one connected component, meaning the disease can get trapped within smaller subpopulations of individuals, limiting its effect.

Discussion
We have proposed a model describing the time evolution of an SIR epidemic spreading through a population of individuals in a multiplex network consisting of two layers: a static network layer representing persistent human connections and a dynamic network layer representing temporary human interactions made outside of a typical household. The model incorporates heterogeneity in the structure, type and duration of connections between individuals, and the number of model equations remains fixed regardless of population size. We designed the multiplex model to afford control of network transitivity (clustering), on the static layer only, by generating the associated network structure using a combination of 2-vertex and 3-vertex cliques, referred to here as static lines and triangles. The dynamic network layer was generated via a single distribution for 2-vertex cliques. Following the EBCM approach (Miller 2014), we obtained expressions for time-evolving quantities of interest, such as the infectious proportion of the population I (t). We have also applied the next-generation matrix method (Diekmann et al. 2009) to compute the basic reproduction number R 0 , a measure of the expected number of infections a typical infectious individual will cause during an epidemic.
Multiplex model equations (1)-(23) were validated, first by testing convergence of epidemic dynamics to predictions made by existing uniplex edge-based compartmental model equations, when either network layer (static or dynamic) was eliminated, and second by comparing full model (with static and dynamic elements) predictions to the dynamics of corresponding statistically correct Gillespie simulations (Gillespie 1976).
The multiplex model's parameter space was explored by varying individual parameters and plotting the resulting epidemic dynamics, and by mapping the outcome on final epidemic size of having various proportions of each edge type when considering different combinations of model parameters β s , β d and η. The basic reproduction number R 0 , found via the leading eigenvalue of next-generation matrix G (24), was validated by demonstrating that continually incrementing infection parameters β s and β d , with all else held constant, led to a rapid increase in final epidemic size as R 0 exceeded its epidemic threshold. Finally, we explored the effect on R 0 and the associated final epidemic size predicted by Eqs. (1)-(23) of altering specific multiplex network attributes governing the rate of rewiring, the extent of clustering and the average degree.
Our unique contribution towards the literature is a model with a combination of static and dynamic network elements, derived by combining the EBCM approach to modelling an SIR epidemic on a static network with tunable clustering (Volz et al. 2011) with the EBCM approach to modelling an SIR epidemic on a dynamic fixeddegree network (Miller et al. 2012), under the framework of a dual-layer multiplex network.
The EBCM approach allows us to model variations in contact structure, contact type and contact duration simultaneously. Modelling such heterogeneities via EBCM provides an opportunity to investigate the effects of heterogeneities observed in realworld networks (Perry-Smith and Shalley 2003;Komurov and White 2007;Vernon and Keeling 2009), alongside consideration of common network attributes such as clustering and degree distributions. EBCM also affords a huge reduction in the number of equations required to track the epidemic, compared with full simulation.
This work progresses the drive to derive population models that capture reasonable levels of complexity and heterogeneity whilst exhibiting a tractable number of equations. By providing a clear and concise 'walkthrough' to deriving and validating our desired model, we hope that future researchers are inspired to build on these results by designing and implementing novel models, modelling approaches and computational algorithms.
The work here extends previous research following the edge-based compartmental modelling approach. Prior EBCM approaches derived model equations describing the SIR epidemic process on wholly static or wholly dynamic uniplex networks. For example, EBCM has been utilised to describe the SIR epidemic on static actualdegree configuration model (CM) networks (Miller et al. 2012), static CM networks with tunable clustering (Volz et al. 2011) and static expected degree mixed Poisson (MP) networks (Miller et al. 2012).
Dynamic uniplex networks have also been considered via the EBCM approach. Namely, CM networks with mean-field social heterogeneity (edges are broken and rewired at a very fast rate, meaning all pairs of individuals contact each other at the same rate, and edge durations are fleeting), dynamic fixed-degree CM networks (edges are rewired, but edge durations are finite), dormant contact CM networks (existing edges are broken and remain dormant for some time, before being re-established), MP networks with mean-field social heterogeneity (fleeting edge duration) and dynamic variable-degree MP networks (finite edge duration) (Miller et al. 2012).
Existing modelling approaches incorporating heterogeneity include the consideration of an epidemic with two 'levels' of mixing between individuals (but no network structure) (Ball et al. 1997), and the later considerations of epidemic processes occurring on structured populations with two levels of mixing , and with two routes of transmission (Zhao et al. 2014). Recently, the EBCM approach was used to derive equations describing an SIR epidemic process with non-sexual and sexual transmission routes, a characteristic of diseases such as Ebola and Zika (Miller 2017).
Other modelling approaches have incorporated dynamicity of connections between individuals (and hence heterogeneity in contact duration) by, e.g. considering an SIR epidemic on a network with intermittent social distancing, where susceptible individuals break links with infectious individuals for some time t b , after which the connection is re-established (Valdez et al. 2013). Another approach considered the effects of constrained rewiring during an SIS epidemic, whereby susceptible individuals cut links to infectious individuals regardless of distance, and rewire to a susceptible individual within a given radius, where the nodes of the network were embedded in Euclidean space (Rattana et al. 2014).
Research considering the large graph limit of an SIR epidemic on a dynamic multilayer network affords heterogeneity in contact type and in contact duration by allowing individual network layers to contain either activating or de-activating edges and by allowing edges in different layers to correspond to different types of contacts (Jacobsen et al. 2016). Although Jacobsen et al. (2016) consider the SIR epidemic spreading on a multiplex network, including providing a dual-layer multiplex example where edge types correspond to community and healthcare contacts, they do not consider any fully static network components.
There are a number of adaptations that can be made to the proposed model. The model considers a heterogeneous contact structure between N individuals. However, the locations of N individuals are not taken into account. Real-world networks occur in space (Barthélemy 2011), and thus, it is important to investigate the effects of considering node locations. In this study, we have chosen to disregard the spatial locations of individuals. A more realistic model of an SIR epidemic spreading on a multiplex network of individuals would be achieved by embedding the locations of each individual into Euclidean space. Even more complex models could consider dynamic node locations or a combination of static and dynamic node locations.
Another potential adaptation is considering weighted network connections. In the proposed model, all connections are considered to be unweighted or equivalently to share equal weight (homogeneity). The model could be adapted by, e.g. making the weight of each connection proportional to the Euclidean distance between the two node locations (given spatial embedding), by imposing a distribution of connection weights or by assigning weights at random. Then, the probability of contracting disease across a connection can be made proportional to the weight of that connection.
In the present model, the population of N individuals is fixed. We do not consider the effect of flux in or out of the population, e.g. by births, deaths or migration events. An important next step is to adapt the model presented here to consider in-and outflow of members of the population or at least to consider whether such in-and outflows significantly influence disease dynamics.
Another model limitation concerns the assumptions made surrounding edge rewire events on the dynamic network layer. Here, we assume that when one partnership ends, a new partnership forms immediately. Thus, given a nonzero degree on the dynamic network layer, individuals remain connected to strangers from the wider population (via dynamic network connections) at all times. In reality, the fleeting connections an individual makes with strangers are temporary, and individuals can remain disconnected from these connections for some time. An improvement to the model could thus be achieved by allowing for gaps to occur between partnerships by implementing the dormant contact approach on the dynamic network layer [e.g. see section 3.3 of Miller et al. (2012), Valdez et al. (2013), Shkarayev et al. (2014) and Tunc et al. (2013)]. The immediate implication of such an approach is a more accurate model in relation to observed human behaviours. However, the dynamics of the epidemic process will be slowed, especially if the duration of the gap (in time) between partnerships is comparable to or longer than the typical time it takes to transmit infection to a partner. Alternative rules for edge dynamicity can also be considered, such as constrained rewiring (Rattana et al. 2014) and edge activation and deletion (Jacobsen et al. 2016;Sélley et al. 2015;Taylor et al. 2012). Other model adaptations include allowing for tunable clustering on all network layers (and thus imposing two edge distributions on each network layer), implementing more complex distributions governing the degrees of each node and biasing initially infectious individuals instead of selecting them at random.
The multiplex model affords tunable clustering on the static network layer by generating its contact structure using a distribution of line stubs and a distribution of triangle corners. However, the configuration model wiring process requires that any two individuals share at most one connection within a single network layer. Double edges can occur across network layers (i.e. when the same edge is present in both network layers), but not within them. This constraint greatly reduces the possibilities for placing triangles suitably into the network, meaning the configuration process is slowed down and the extent of clustering that can be achieved is reduced. Greater control over clustering could be achieved by adapting the model to allow for overlapping triangles (and either allowing double edges to occur in single network layers, or amalgamating any double edges that occur into single edges, or doubly weighted edges).
Other than making adaptations to the proposed model, there are a number of tests and analyses which are beyond the scope of this work. Firstly, a comprehensive exploration of the entire parameter space would elucidate the behavioural 'envelope' of the model and uncover any parameter regions where the model poorly predicts the SIR epidemic process, compared with simulation. A more thorough understanding of the impact of degree and degree heterogeneity on the relationship between parameters and system behaviour will require consideration of additional edge distributions with various levels of heterogeneity and average degrees. Secondly, the model's utility can be investigated by using real-world data from historical epidemics or similar processes, e.g. livestock herd contact tracing data or Twitter data tracking the prevalence of a hashtag over time. Using real data, model parameters could be estimated using Bayesian estimation techniques and the resulting model predictions compared with prior knowledge of what occurred. The basic reproduction number R 0 can be tested in the same way.
This work considers an SIR compartmental model under the guise of a disease spreading through a networked population. Thought must be given to what other realworld processes can be well described by the SIR compartmental model, such as opinion formation, rumour spreading or uptake of fashion trends. Further, a two-layer multiplex like the proposed model could be used to investigate the dynamics of two interacting SIR-type processes, such as a physical disease spreading process occurring on one network layer in combination with a disease awareness process occurring on the opposing network layer, using similar approaches to those of Funk et al. (2009) and Li et al. (2015).
Future research can build on these observations by considering similar modelling approaches that account for compartmental models other than the SIR type. For example, the SIS model (describing infections that do not confer lasting immunity, such as the common cold) and the SEIR model (describing infections with incubation periods, where individuals have contracted a disease but are not yet infectious and hence are in the 'exposed' disease state) are not considered here. Modelling an SEIR infection may require simple adaptation of the existing EBCM approach. However, consideration of an SIS-type epidemic process requires an altogether new modelling approach.
A key assumption of the present approach is the consideration of all neighbours of the test node u as being independent. Attempting to impose this assumption would prevent modelling of SIS dynamics, a consequence which is discussed in Miller and Kiss (2014) and Miller et al. (2012).
Experiments that can be performed to improve and inform future modelling approaches include: quantifying the levels of heterogeneity in existing populations, including behavioural and structural heterogeneity, gaining a deeper understanding of the biological processes underlying disease spreading processes, improving on existing algorithmic and analytic approaches and fostering closer relations between modellers and practitioners, in order to maximise the benefits arising from research.