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 [1].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 [2,3,4,5], and further that an epidemic can spread on a multiplex network even if the individual layers are well below their respective epidemic thresholds [6].A global cascades model generalised for multiplex networks was used to show that multiplexes are more vulnerable to global cascades than single layer networks [7].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 [8].
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 [9].Network quality functions able to detect community structure in multiplex networks have been developed [10].Further, results such as the large graph limit of an 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 [11].
In network models, increased clustering is generally considered to slow an epidemic by increasing the epidemic threshold [12].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 [13].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 [3].However, the authors also discovered a critical threshold for the average degree, above which clustering was shown to facilitate global cascades [3].A uniplex network study found that simultaneously increasing clustering and the variance of the degree distribution led to an increase in final epidemic size [14].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 [15].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 [16] 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 [17], random link activation and deletion [18,19,20], and temporary link deactivation [21,22].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 fixed-and dynamic variable-degree models, where edges are swapped at a given rate, or edges are broken and created at given rates, respectively [23,24].
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 [25,26,27,24], 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, contactduration, and contact-structure.We derive the associated basic reproduction number R 0 , following the next generation matrix approach [28].We describe the implementation of the EBCM model and of statistically-correct Gillespie simulations of the epidemic process [29].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 the 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 , the expected number of static line stubs that belong to a randomly selected individual, M , the expected number of static triangle corners that belong to a randomly selected individual, and M , 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 [24].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 , 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 have 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 is shown in , and using the initial condition φ R (0 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 [30].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− The proportion of individuals in the network that are a member of s static line stubs, t triangle corners and d dynamic line stubs g(x, y, z) Probability generating function for the numbers of static lines, triangles and dynamic lines of which an individual is a member θ 2 (t) 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 φ S (t),φ 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 A(t) 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 B(t) 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 ψ S (t),ψ I (t),ψ R (t) 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 Table 1: Definitions for model variables and parameters.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.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 ).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 4 (recall that u cannot transmit infection), where s, t and d are realisations of the excess degree distribution.We calculate φ S by multiplying the probability that a random neighbour across a static line has (s, t, d) neighbours, with the probability the random neighbour is susceptible, summing over all possible values of (s, t, d), and dividing by M = 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 equations ( 2)-( 4), and noting from Fig Considering θ 3 Since θ 3 denotes the probability that neither of the other nodes in a triangle have 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 have 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 satisfies 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 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 = −2Aφ SS , which 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: 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.

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 have 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 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.The flux between π S , π I and π R , the probabilities that a randomly chosen dynamic stub belongs to a susceptible, infected or recovered node, respectively.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 re-wiring, 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 s,t,d / M , where M = g (z) (1, 1, 1) 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 M , 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 The equation for π S can be condensed using the 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 dynamic-edge 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 , 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 have ζ = ηθ 2 4 − ηζ.We have an expression for ζ, so the derivative of ψ S can be found via the chain rule: with simplifications achieved by utilising π 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 duallayer 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 ij 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 [28], the value of R 0 is found via the leading eigenvalue of the matrix G, or equivalently, the eigenvalue with greatest magnitude.
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 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 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 2βs 2βs+γ 2βs 2βs+2γ = 2βs 2βs+γ βs βs+γ .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 Then we have where k t denotes the expected static triangle corner degree.Finally, we have , where t denotes the expected excess static triangle corner degree.
We conclude by deriving elements from the third column of G, starting with non-diagonal 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 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 r = 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 = . The basic reproduction number R 0 is the eigenvalue of the 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 edge-based compartmental model type, and of a user specified end time for the computation.
Solutions to equations ( 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 equations ( 1)- (23), found in the MATLAB environment, Gillespie simulations [29] 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 are 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 ensures that the total number of static line stubs is even, the total number of dynamic line stubs is even, and that the total number of static triangle corners is a multiple of three.
Firstly, the static network layer is generated using vectors containing the number of static line stubs and triangle corners each individual partakes in.Pairs of static line stubs and triples of static triangle corners are selected at random.Provided potential static lines and triangles do 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 are added to the static adjacency matrix.The unmatched static line stubs and static triangle corners lists are updated, and the process continues until all static line stubs and triangle corners are successfully matched.
Secondly, the initial structure of the dynamic network layer is generated using the vector storing the number of dynamic line stubs each individual partakes in.Pairs of dynamic line stubs are selected at random.Provided a potential dynamic edge does not generate a self-loop or a doubleedge within the dynamic network layer, it is added to the dynamic adjacency matrix.Successfully paired dynamic stubs are removed from the unmatched stubs list, and the process continues until all dynamic line stubs are successfully matched.
The nature of this configuration model approach means the wiring processes for the static and dynamic network layers may have to be restarted multiple times in order to achieve final network structures.Once all static line stubs, static triangle corners and dynamic line stubs have been wired up, the configuration process is complete.Although the script prevents double-edges from occurring within each network layer, it is 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 can be implemented.In each Gillespie simulation, ρN initially infectious individuals are selected at random from the population.At each time step, a vector of length (N + 1) describes 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 follow an exponential distribution with scale parameter 1 R , where R := the sum of the rates vector at the current time step.Each event occurring is either an infection, a recovery or an edge swap.Uniformly distributed random numbers are generated at each time step to determine the next event to occur.When an edge swap event occurs, the script selects two dynamic edges at random, ensuring that all four nodes involved in these edges are unique.The script also ensures that the proposed new dynamic edges do not already exist within the dynamic network layer.Given these conditions, an edge swap occurs and the Gillespie process continues.The process terminates once the user specified end time is reached.Matlab code written to implement the above routine will be provided upon acceptance of the manuscript.

Results
We have followed the EBCM approach to derive equations ( 1)-( 23), which describe an SIR epidemic spreading on 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 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.We have also applied the next generation matrix method [28] 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.
In what follows, we assess the validity of equations ( 1)-( 23) and of the basic reproduction number R 0 , obtained via the 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 equations ( 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 equations ( 1)-( 23) in predicting the epidemic process for a number of multiplex network configurations.By solely considering the predictions of equations ( 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 equations ( 1)-( 23) when altering the rate of rewiring, the extent of clustering, and the average degree in the multiplex model.

Model convergence to existing uniplex model equations
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 [14].
By comparing predictions made by uniplex model equations in [14] with those of multiplex model equations ( 1)-( 23) when dynamic network elements are 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 [14] 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 [24].
Excellent agreement was observed between predictions made by equations ( 1)-( 23) when static network elements are close to zero, output from the dynamic fixed-degree derivation in [24] and Gillespie simulations describing the SIR epidemic and edge rewire processes occurring simultaneously on equivalent multiplex networks, for a number of setups 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 equations ( 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 equations ( 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 generate a negative binomial distribution with smaller average degree and a reduction in variance, slowing the epidemic's spread.Larger values of r 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 ) is used to plot dynamics predicted by multiplex model equations ( 1)-( 23) (thick black line).In each panel, a single parameter is varied and the resultant predictions are plotted in various colours, indicated by individual panel legends.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  2 the difference is added to each of the values p t and p d .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 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
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 equations ( 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
The 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 sub-threshold (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 equations ( 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 sub-threshold values, while 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 equations ( 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.Again, the final epidemic size begins to increase more quickly after some critical average degree has been reached, around k = 14.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 sub-populations 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 dual-layer static-dynamic multiplex network.The model incorporates heterogeneity in the structure, type and duration of connections between individuals.Following the EBCM approach [27], we obtained expressions for time-evolving quantities of interest, such as the infectious proportion of the population I(t).An estimate of the associated basic reproduction number R 0 was derived, utilising the next generation matrix method [28].
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 statisticallycorrect Gillespie simulations [29].
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 the 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 equations ( 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 [14] with the EBCM approach to modelling an SIR epidemic on a dynamic fixed-degree network [24], 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 real-world networks [31,32,33], 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 actual-degree configuration model (CM) networks [24], static  CM networks with tunable clustering [14], and static expected degree mixed Poisson (MP) networks [24].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) [24].
Existing modelling approaches incorporating heterogeneity include the consideration of an epidemic with two 'levels' of mixing between individuals (but no network structure) [34], and the later considerations of epidemic processes occurring on structured populations with two levels of mixing [35], and with two routes of transmission [6].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 [36].
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 [37].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 [17].
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 [11].Although Jacobsen et al. [11] 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 [38] 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.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, biasing initially infectious individuals instead of selecting them at random, and considering alternative rules for edge dynamicity, such as constrained rewiring [17], edge activation and deletion [11,19,18], or the dormant contact approach [37,22,21].
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 real-world 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 [39] and [40].
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 [41] and [24].
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.

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 can occur prior to time zero.By integrating the relation dφR dt = γ βs d(1−θ2) dt

Figure 1 :
Figure 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 ).

Figure 2 :
Figure 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.

Figure 3 :
Figure 3: 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.

Figure 4 :
Figure 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.
Fig 3 we have ψS = ηθ 4 π S − ηψ S − Bψ S , so we calculate the flux between compartments ψ S and ψ I using the rate

1 −r + d β d β d +γ 1 1
ηβ d (η+γ)(β d +γ) .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 β d +γ 1 1−r .Thus we find G dd = r −r , where d is the expected excess dynamic line stub degree.

Figure 5 :
Figure 5: Multiplex model convergence -no dynamic layer, with simulation.The time evolution of infection prevalence for the original EBCM of an SIR epidemic on a static uniplex network (solid black line), for the proposed EBCM of an SIR epidemic on a dual-layer multiplex with the dynamic network layer being close to zero (thick dashed red line), and for 10 Gillespie simulations of the SIR epidemic on a single network of size N = 5000 (solid blue lines).In all panels γ = 1, ρ = 0.05, and p = 0.5 and r = 10 generate a negative binomial distribution for pairs of edge stubs.For the original static derivation (solid black line) p s = 0.5 = p t , describing the proportion of edge-pairs that are split into two single lines or remain as a triangle corner, respectively.For the multiplex derivation (thick dashed red line) β s = β d , η = 0.01, and p s = 0.4999999, p t = 0.5 and hence p d = 10 −7 describe the proportion of edge-pairs that become two static lines, a static triangle corner, or two dynamic edges respectively.(A) β ′ s = 1, C = 0.02677, (B) β ′ s = 0.5, C = 0.02670, (C) β ′ s = 0.25, C = 0.02658, (D) β ′ s = 0.125, C = 0.02685, where C denotes the global clustering coefficient of each static network layer generated for simulation.

Figure 6 :
Figure 6: Multiplex model convergence -no static layer, with simulation.The time evolution of infection prevalence for the original EBCM of an SIR epidemic on a dynamic uniplex network with conserved degrees and edge re-wiring (solid black line), for the proposed multiplex EBCM of an SIR epidemic with the static network layer being close to zero (thick dashed red line), and for 10 Gillespie simulations of the process on a single network of size N = 5000 (solid blue lines).In all panels γ = 1, ρ = 0.05, and p = 0.5 and r = 10 generate a negative binomial distribution for pairs of edge stubs.For the original conserved-degree derivation (solid black line) p d = 1, indicating that all edge-pairs become two disjoint dynamic edges.For the multiplex derivation (thick dashed red line), η = 0.01 and p s = p t = 10 −7 and p d = 0.9999998 describe the proportion of edge-pairs that become two static lines, single triangle corners, or two dynamic edges respectively.(A) β ′ s = 1, C = 0.004944, (B) β ′ s = 0.5, C = 0.005285, (C) β ′ s = 0.25, C = 0.005344, (D) β ′ s = 0.125, C = 0.005127, where C denotes the global clustering coefficient of each dynamic network layer generated for simulation, at time zero.

Figure 7 :
Figure 7: Multiplex model prediction vs. simulation -varying clustering and average degree.Plotting the dynamics of the proportion of infected individuals over time.Each panel contains 25 Gillespie simulations on a single multiplex network comprised of N = 1000 individuals (blue lines) and the associated EBCM prediction (black line).All networks are generated using a negative binomial distribution for pairs of edge stubs with parameters p = 0.5 and various values for r.Networks in column 1 (counting from left to right) have average degree 10 (achieved via r = 5), networks in column 2 have average degree 20 (achieved via r = 10), and networks in column 3 have average degree 30 (achieved via r = 15).Networks in row 1 (counting from top to bottom) have minimised clustering via values p s = 0.99999998 and p t = 10 −8 .Networks in row 2 have the values p s = 0.49999999 = p t .Networks in row 3 have maximised clustering via the values p s = 10 −8 and p t = 0.99999998.Counting panels from left to right and top to bottom, starting with the upper-left panel, static networks have the following clustering coefficients: C = 0.0161, C = 0.0267, C = 0.0370, C = 0.0535, C = 0.0473, C = 0.0493, C = 0.0898, C = 0.0662, C = 0.0629.In all panels, tmax = 10, ρ = 0.05, β s = β d = 0.25, γ = 1, η = 0.01.

Figure 8 :
Figure 8: Multiplex model prediction vs. simulation -varying infection parameters β s and β d .Plotting the dynamics of the proportion of infected individuals over time.Each panel contains 100 Gillespie simulations (10 simulations on 10 multiplex networks comprised of N = 5000 individuals) (blue lines) and the associated EBCM prediction (black line).All multiplex networks follow a negative binomial distribution for pairs of edge stubs with parameters p = 0.5 and r = 10, which were split into three edge types via p s = 0.3 = p t and thus p d = 0.4.In all panels tmax = 10, ρ = 0.05, γ = 1, η = 0.01.Across the panels, different values for β s and β d have been used in the range [0.125, 0.25, 0.5], indicated by individual column and row headings.

Figure 9 :
Figure 9: Multiplex model predictions.Plotting the dynamics of the proportion of infected individuals over time, for a number of different parameter sets.In all panels, a baseline parameter set (p = 0.5, r = 10, p s = 0.3 = p t , p d = 0.4, β s = 0.05,β d = 0.2, γ = 1, η = 0.01 = ρ, tmax = 10 ⇒ R 0 = 1.076) is used to plot dynamics predicted by multiplex model equations (1)-(23) (thick black line).In each panel, a single parameter is varied and the resultant predictions are plotted in various colours, indicated by individual panel legends.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 subtract1  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  2 the difference is added to each of the values p t and p d .

Figure 11 :
Figure 11: Validation of the basic reproduction number R 0 .Plotting values of the basic reproduction number R 0 (x-axis), found via the leading eigenvalue of the matrix (24), against the associated final epidemic sizes (y-axis) predicted by multiplex equations (1)-(23) (red line) and recorded by single statistically-correct Gillespie simulations (blue circles).Static and dynamic line stubs follow Binomial distributions with parameters n = 20 and p = 0.5.The distribution of triangle corners follows a Binomial distribution with parameters n = 1 and p = 0.001 to minimise clustering.Fixed parameters were γ = 1, ρ = 0.001, η = 0.01, tmax = 10, N = 1000.In each setup β s = β d .100 transmission rates were tested, from β s = β d = 0.01 up to β s = β d = 0.3, in equal-sized increments.In Gillespie simulations where R 0 > 1, if the number of infected individuals did not reach 10 times the initial number of infectives, all data was discarded and the Gillespie script restarted from initial conditions at time zero.