Dynamics of epidemic models from cavity master equations
aa r X i v : . [ phy s i c s . s o c - ph ] J un Dynamics of epidemic models from cavity master equations
Ernesto Ortega, David Machado, Alejandro Lage-Castellanos
We apply the cavity master equation (CME) approach to epidemics models. We explore mostlythe susceptible-infectious-susceptible (SIS) model, which can be readily treated with the CME as atwo-state. We show that this approach is more accurate than individual based and pair based meanfield methods, and a previously published dynamic message passing scheme. We explore averagecase predictions and extend the cavity master equation to SIR and SIRS models.
I. INTRODUCTION
Since the seminal works introducing susceptible-infectious-recovered compartment models (SIR) of Kermack andMcKendric [1], epidemics modeling has grown fast as a field. The approach has changed with time, from dynamicalsystems or population dynamics towards a more stratified approaches as patchy, mobility based, age-structured orcontact matrices compartment models.The current context of a global COVID-19 outbreak and the perspective of a coexistence with an endemic virusrequires a test-trace-isolate epidemiological system to keep the outbreak controlled. Much attention is now put onagent based models that could improve the efficacy of the testing strategy. Assuming that new technologies canprovide reliable contact data between humans, the likelihood of people being infected needs to be estimated either bynumerical simulations or some statistical modeling. To this end, it is suitable to count with fast algorithms that canaccurately predict probabilities of infection for agents in networks.There are a wide variety of such algorithms. The classical approach to the forecasting of epidemics on networks is anaveraging of the master equation of the process complemented by a factorization assumption at some level. This yieldsa hierarchy of ever more complex but more accurate differential equations for expected values and correlations[2, 3].Most of the time only the first two levels, known as individual based mean field and pair based mean field, are used.More recently, ideas from discrete optimization algorithms have sneaked into the inference of SIR kind of modelsin the shape of Dynamical Message Passing [4, 5] or Belief Propagation (BP) [6, 7]. The main difference with respectto the previous approach is the appearance of conditional probabilities, rather than correlations, to be integrated intime. It has been used with success in the reconstruction of epidemics on graphs and it is currently being tried in thetask of risk assessment for COVID-19.It is known that BP’s fixed point is connected to the cavity method from statistical mechanics [8]. An extensionof the latter to continuous time Markov Chain processes within discrete spin systems have been recently achievedthrough the derivation of a set of differential equations for cavity conditional probabilities: the cavity master equation(CME) [9].In this article we explore CME’s application to SIS and SIR-like models in graphs (section III). When comparedwith Monte Carlo simulations of the epidemic, it is shown to outperform three other approximations. In some simplecases we draw analytical results for the steady state and critical reproduction number (section IV). We also explorethe average case dynamics of these equations, in particular in random regular and Erdos-Renyi graphs (sec V). Weconclude by briefly showing the application of CME to SIR and SIRS models in section VI.
II. EPIDEMICS ON NETWORKS
In what follows we focus on continuous time compartment epidemic models on networks. We assume a fixed networkof contacts to be given G = ( V, E ) with a set of vertices V = { , . . . , N } and a set of edges E . An edge ( i, j ) ispresent if nodes i and j are neighbors in the network, meaning there is a possibility of transmission of diseases betweenboth nodes.The standard susceptible-infectious-susceptible model (SIS) considers the nodes to be in either of two compartments(states) X i = 0 ≡ susceptible, or X i = 1 ≡ infectious and is the simplest standard for recurrent transmissible diseases.The epidemic is thus a continuous time stochastic process with only two admitted transitions occurring at • a rate β , at which a link ( i, j ) can transmit the state 1 from node i to j ; • a rate µ at which state 1 decays to state 0 on any infectious node.as represented in diagram 1. An analytical description of this stochastic process is given by the master equation forthe evolution in time of the probability over the whole configuration space P ( X i , . . . , X N , t ) [3]. However concise andexact, the integration in time of such equation is generally impractical given the size 2 N of the configuration space. Fig. 1.
Allowed transitions in SIS compartment model on networks.Attempts to reduce the complexity start from factorizing the single master equation into many equations for eachnode marginals P i ( X i , t ). Given that the X ’s are two state variables, P i ( X i ) is parameterized by the mean valueE[ X i ] = P i ( X i = 1). This results in an equation that is still exact [3]dE [ X i ( t )]d t = E − µX i ( t ) + (1 − X i ( t )) β N X j =1 a ij X j ( t ) (1)where a ij = 1 if ( i, j ) ∈ E and is zero if nodes i and j are not neighbors in the graph.However, the expectation value on the right hand side acts over products of variables X i ( t ) X j ( t ), which requiresa differential equation for the evolution of the correlations. Not surprisingly, the two point correlations functionsdepend on three point correlations and so on and so forth.The simplest closure of equation eq. (1) is the individual-based mean field (IBMF) in which independence is assumedE[ X i ( t ) X j ( t )] = E[ X i ( t )] E[ X j ( t )] ≡ ρ i ( t ) ρ j ( t ) and therefore eqs. (1) are now a closed set of non linear differentialequations: dρ i ( t ) dt = − ρ i ( t ) + λ [1 − ρ i ( t )] N X j =1 a ij ρ j ( t ) (2)The second simplest closure is that known as pair-based mean field (PBMF), in which two point correlations aretreated analytically [10]:dE [ X i X j ] dt = − µ E[ X i X j ] + β P Nk =1 a ik E[ X j X k ] + β P Nk =1 a jk E[ X i X k ] − β P Nk =1 ( a ik + a jk )E[ X i X j X k ] (3)but a factorization is assumed for higher correlations. Different approaches have been used to approximate E[ X i X j X k ]in terms of smaller correlations. We will compare to that in which E[ X i X j X k ] ≡ E[ X i X j ]E[ X k ] ≡ ρ ij ( t ) ρ k ( t )[10]d ρ i ( t )d t = − µρ i ( t ) + β X j ∈ ∂i [ ρ j ( t ) − ρ ij ( t )] (4)d ρ ij ( t )d t = − µρ ij ( t ) + β [ ρ i ( t ) + ρ j ( t ) − ρ ij ( t )] ++ β X k ∈ ∂i \ j [ ρ jk ( t ) − ρ ij ( t ) ρ k ( t )] + β X k ∈ ∂j \ i [ ρ ik ( t ) − ρ ij ( t ) ρ k ( t )] (5)In both approaches, IBMF and PBMF, the expected values evolving in time are intended to be expectations overdifferent stochastic stories of the whole epidemic process. Therefore they are to be compared with averages over manyMonte Carlo simulations of such process.A slightly different approach to modeling epidemics on graphs come from message-passing inspired methods. Thedynamic message passing [4, 5] involves a set of probabilities and a set of conditional probabilities, rather thancorrelations [5]: d p i d t = − µp i + β (1 − p i ) X k p ki (6)d p ij d t = − µp ij + (1 − p j ) β X k ∈ ∂i \ j p ki (7)where p i = E[ X i ] just as before, but p ij ≡ P ( X i = 1 | X j = 0) is a conditional probability resembling the kind ofcavity fields or messages that commonly appear in the cavity method and the belief propagation algorithm for discreteoptimization.The main results of this article are drawn from the cavity master equation [9]. They are quite similar to those ofrDMP and help formalizing this rather empirical approach by deriving it from a more solid mathematical setting. Indoing so we not only correct one term of the rDMP equations, but also underline the approximations involved andtherefore shed light on possible improvements. III. CAVITY MASTER EQUATIONS FOR SIS EPIDEMICS
In order to connect with its first presentation in [9], we start by considering the general continuous time dynamicsof a system σ = { σ , . . . , σ N } of N bimodal variables σ i ∈ {± } interacting with their nearest neighbors in somegiven topology. In the very generic Markovian case, the dynamic is fully defined by the rate function r i ( σ ) at whichvariables flip their states from σ i → − σ i . The distribution P ( σ , t ) in the configuration space of this stochastic processis ruled by the joint master equation: dP ( σ ) dt = − N X i =1 h r i ( σ ) P ( σ ) − r i ( F i ( σ )) P ( F i ( σ )) i , (8)where F i represents the flip operator on variable i , i.e. F i ( σ ) = { σ , . . . , σ i − , − σ i , σ i +1 , . . . , σ N } .Although exact, the previous equation is useless already for middle size systems, since it actually represents a set of2 N coupled differential equations that take exponential time to enumerate, let alone to integrate. However this is thecorrect starting point for approximations, as is usually done to obtain mean field approximations (2,4,5) in epidemicsmodels [2, 3, 11].In [9] this master equation is recast into an equilibrium problem by extending the configuration space to considerthe continuous trajectory of each variable in time X = { X , . . . , X N } where X i = { σ i ( t ) : ∀ t ∈ [ o,T ] } . Althoughat glance it seems untreatable the infinite dimensional space for the functions X i ( t ), the discrete values of σ i ( t )allows for a codification of the functions in a numerable number of transitions times { T { i } , T { i } , . . . } such that σ i ( T { i } k ) = − σ i ( T { i } k + d t ). The resulting Random Point Process is treated with standard techniques in equilibriumstatistical mechanics to write down a closed set of cavity master equations as:d P ( σ i )d t = − X σ ∂i h r i ( σ i , σ ∂i ) (cid:2) Y k ∈ ∂i p ( σ k | σ i ) (cid:3) P ( σ i ) − r i ( − σ i , σ ∂i ) (cid:2) Y k ∈ ∂i p ( σ k | − σ i ) (cid:3) P ( − σ i ) i (9)d p ( σ i | σ j )d t = − X σ ∂i \ j " r i [ σ i , σ ∂i ] h Y k ∈ ∂i \ j p ( σ k | σ i ) i p ( σ i | σ j ) − r i [ − σ i , σ ∂i ] h Y k ∈ ∂i \ j p ( σ k | − σ i ) i p ( − σ i | σ j ) (10)We have lighten the notation by not putting the i and ij dependence of the distributions, understanding that theyassume the index of the variables they depend on (as P ( σ i ) ≡ P i ( σ i )). Both probabilities are intended in the sense“over the ensemble of dynamic evolutions up to time t ”, starting from the same initial conditions.This is a substantial improvement over the original master equation since we are dealing now with N functions P i ( σ i , t ) representing the distribution of variables σ i , and with N ∗ h k i functions p ( i,j ) ( σ i | σ j , t ) that represent theprobability of finding variables in a given state, conditioned to the state of one of its neighbors ( h k i is the averagedegree of the nodes in the interactions network). In the worst case of a fully connected system, we still would have O ( N ) equations, that can be numerically integrated even for relatively large systems, compared to what we can dowith eq. (8).We can readily translate the cavity master equations (9,10) to the case of S.I.S epidemic model by identifying ourtwo states as S → σ i = − I → σ i = 1. After complementary ( P i ( I ) + P i ( S ) = 1), it is enough to track theinfection probability in each node P i ( I ). The first term in equation (9) largely simplifies due to the fact that therecovery rate r i ( σ i = 1 , σ ∂i ) = µ is independent of the neighbors state, resulting in − µP ( I ). Considering also thatthe transmission rates are additive r i ( S, σ ∂i ) = P k ′ r i ( S, σ k ′ ) and that r i ( S, S ) = 0 and r i ( S, I ) = β we get to thecavity master equations for the S.I.S model asd p i d t = − µp i + β (1 − p i ) X k p ki (11)d p ij d t = − µp ij + (1 − p ij ) β X k ∈ ∂i \ j p ki (12) P r ob -I n f- t MC simCMErDMPIBMFPBMF L - E rr o r t CMErDMPIBMFPBMF Fig. 2. Left:
Probability of node 29 to be infected as a function of time, discounting at each time the situations inwhich the epidemics disappear. The epidemic outbreak was in node 1, with β = 0 . µ = 0 .
05. The inset is acloser look in the stationary state.
Right:
L1 distance between the four approximations IBMF, PBMF, rDMP,CME with respect to MC simulations.where we simplified notation further by making P ( σ i = 1) ≡ p i and P ( σ i = 1 | σ j = − ≡ p ij . Appendix A containsa more detailed derivation of these equations for SIS model.Equations (11) and (12) are almost identical to those obtained for rDMP (6) and (7), except for the very lastterm where (1 − p j ) is replaced by (1 − p ij ). Not surprisingly, the results from both approximations might not betoo different in certain cases. For instance, by reproducing the very same case of [5] of an SIS epidemic outbreak inthe Zacharia’s karate club network, starting at node 1 both approaches are quite similar as shown in figure 7 (left)(actually, in that figure is almost impossible to distinguish one from the other and CME seems to be missing). Infigure 7 (right) the L1 distance between the average from many Monte Carlo simulations and the predictions madeby all four methods CME, rDMP, IBMF and PBMF shows features that will repeat in other benchmarks: • that both rDMP and CME get the general qualitative behavior well, • with a rather faster outbreak expansion in the transient and • with CME better fitting the stationary state.Individual based mean field tends to be the fastest growing prediction, while pair based mean field is the closest toCME, yet not equal to it.When computing the Monte Carlo averages we have neglected the simulations in which the epidemic is randomlywiped out in the first few iterations. None of these methods can take this fluctuations into account, since they are allmean-field approaches.Differences are more evident in the case of random regular graphs with degree k = 3 in figure 3. We represent theaverage epidemic size, i.e. the fraction of infected nodes respect to number N = 1000 of nodes in the graph at everytime step. We started from a fraction α = 0 . IV. ENDEMIC (STEADY) STATE
We can easily compute the epidemic size in the endemic state for these approximations in the case of randomregular graphs. Since every node has the same amount of neighbors k , the equations are the same for every node, andwe can assume that in the steady state the topology is averaged out, and all the probabilities are the same, regardlessthe node indexes.Working explicitly for the CME approximation, we set d p i d t and d p ij d t to 0 in equations (11) and (12): µp i = β (1 − p i ) X k p ki ≡ β (1 − p i ) kp ij µp ij = (1 − p ij ) β X k ∈ ∂i \ j p ki ≡ (1 − p ij ) β ( k −
1) (13)where now the indices i and i, j are generic. Solving this system of equations we get for the two variables p i and p ij : p i = 11 + µkβ (1 − µ ( k − β ) p ij = 1 − µ ( k − β (14) P r ob - i n f t MC simCMErDMPIBMFPBMF 0 0.1 0.2 0.3 0.4 0.5 0.6 0 50 100 150 200 250 300 350 400 450 500 L - E rr o r t CMErDMPIBMFPBMF P r ob - i n f t Ro = 0.5 0.1 1 10 100 t Ro = 0.7 0.1 1 10 100 t Ro = 0.9 0.1 1 10 100 t Ro = 1.1MC simCMErDMPIBMFPBMF
Fig. 3. Random regular graph with connectivity k = 3: Top figures - Percent of the population in the infectedstate (left) and the L distance (right), in a graph of 100 nodes with parameters β = 0 .
11 and µ = 0 . Bottom figures : Percent of the population in the infected state for differentvalues Ro in a graph of 1000 nodes, starting with half of the population in the infected state.This can be already tested against simulations of random regular graphs. Furthermore we can obtain the criticalreproduction number R = β/µ above which there is an endemic outbreak (sustained in time epidemic) by solving p i ( k ) = 0, resulting in R ∗ = 1 / ( k − p i at the steady stateas a function of R . When compared to Monte Carlo results, mean field approximations (IBMF and PBMF) give anoverestimation of the epidemic size at small R , and rDMP an underestimation at large R . Meanwhile, CME seemsto be a better prediction mixing the large R behavior of IBMF and PBMF with the small R behavior of rDMP. V. AVERAGE CASE
The systems of equations defining IBMF, PBMF, rDMP and CME could be large and delicate to solve, althougha simple numerical integration normally works. However, in many cases we are interested in general predictions forcertain families of graphs or topologies. In this section we derive an average version of the CME approximation tocharacterize SIS epidemics on random graphs.The simplest description of a graph ensemble is given by the distribution of degrees of its nodes. Since the CMEequations depend on the information coming from the neighbors in the network, it is expected that nodes moreconnected will have a different behavior than those less connected. We therefore attempt to reduce the number ofequations in our system by characterizing all nodes with the same degree by a couple of average parameters p γ = 1 N γ X i : d i = γ +1 p i p γ → = 1 M γ X i : d i = γ +1 X j ∈ ∂i p γij (15) Approximation Endemic state (equilibrium) R ∗ IBMF p i = 1 − ρkλ R ∗ = k PBMF p i = ( R ) k + R − ( k − )( R ) k +2 R R ∗ = p f kk − − k rDMP p i = 1 − µ ( k − β p ij = ( k − β − µkβ R ∗ = k − CME p i = µkβ (1 − µ ( k − β ) p ij = 1 − µ ( k − β R ∗ = k − TABLE I.
Critical reproduction number under four approximations: IBMF, PBMF, rDMP and CME P r ob - s t ea dy - s t a t e R MC simCMEPBMFrDMPIBMF 0 0.7 0.3 1 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 3 4 5 6 7 8 9 10 R o* k PBMFCME rDMPIBMF Fig. 4. Left : Comparison between Monte Carlo simulations and the predicted endemic states of each method as afunction of R . Rigth : Comparison between the epidemic threshold predicted by the methods as a function ofgraphs connectivity.In both cases the normalization factors count the number of terms in the sums: N γ is the number of nodes withconnectivity γ + 1 while M γ is the number of graph’s edges that contain one of these N γ nodes. Averaging theequations (11) and (12), and after some simplifications, we get the average CME:˙ p γ = − µp γ + β ( γ + 1)(1 − p γ ) X γ ′ g link ( γ ′ ) p γ ′ → (16)˙ p γ → = − µp γ → + βγ (1 − p γ → ) X γ ′ g link ( γ ′ ) p γ ′ → (17)where g ( γ ) is a contact degrees distribution. A node extracted randomly from the set of nodes V has degree distribution P ( γ ). However, in order to average the CME equations we rather need to know the degree distribution g ( γ ) of nodesthat are sampled by randomly picking up an edge ( i, j ) ∈ E . Both distributions are related, assuming there is nofurther correlations in the graph ensemble, by: g link ( γ ) = ( γ + 1) P ( γ ) P γ ( γ + 1) P ( γ ) . (18)A very simple case is that of regular graphs, where degree distribution is deltaic P ( γ ) = g ( γ ) = δ k − ,γ . Equations E p i d e m i c - S i ze t MC simCMEav-CME Fig. 5.
Epidemic size as function of time. The epidemic outbreak was in 5 nodes of the graph of 100 nodes andconnectivity 5 for each node, with β = 0 . µ = 0 . p k − ( t ) ≡ p ( t ) and p k − → ( t ) ≡ p → ( t ):˙˜ p = − µ ˜ p + β k ˜ p → (1 − p ) ˙˜ p → = − µ ˜ p → + β ( k −
1) ˜ p → (1 − p → ) (19)whose numerical integration can be compared with Monte Carlo simulations of epidemics in graphs with the samevertex degree k , and with the corresponding integration of the single instance CME equations (11) and (12). Figure5 shows that the steady state is well predicted, while the transient is not, even in comparison with CME itself. Thisis a natural consequence of the lost of the spatial structure in the average case. A. General graph ensembles
Equations (16) and (17) are already a reduction of N differential equations to K average equations, where K is themaximum degree in the graph. However, K itself could be large. We can further simplify by averaging now over thenodes degree and reducing to only two parameters ˜ p → = P γ g link ( γ ) p γ → and ˜ p = P γ P ( γ ) p γ :˙˜ p = − µ ˜ p + β ˜ p → X γ ( γ + 1) P ( γ )(1 − p γ ) (20)˙˜ p → = − µ ˜ p → + β ˜ p → X γ γg link ( γ )(1 − p γ → ) (21)These equations are still not closed, since the right hand sides still depend on the degree based parameters. Theproduct by ( γ + 1) and γ inside the sums in the right hand sides do not allow for a direct connection with thedefinitions of ˜ p and ˜ p → . We will show, however, that in the case of Erdos-Renyi graphs such connection can beobtained, though through some approximations and ansatzs. B. Closure on Erdos Renyi graphs
For an Erdos-Renyi graph node degrees are Poisson-distributed: P ( γ ) = e − κ κ γ +1 ( γ +1)! and g link ( γ ) = P ( γ − κ is the average degree. We can connect the terms inside the sums in (20) and (21) with the derivatives of ˜ p and ˜ p → with respect to the parameter κ by noting that ∂ ˜ p∂κ = − ˜ p + κ P γ ( γ + 1) P ( γ ) p γ (22) ∂ ˜ p → ∂κ = − ˜ p → + κ P γ γ g link ( γ ) p γ → (23) f t MC sim
CMEav-CMER = 40.0R = 4.0R = 2.0R = 1.3 Fig. 6.
Comparison between average-CME and CME for Erdos Renyi graphs with mean connectivity c = 3,several values of µ and β = 0 .
4. The figure shows the single site probability of infection as a function of time. As theratio β/µ decreases, the steady state has less infection probability.As it is shown in appendix B, by substitution in (20) and (21) we get:˙˜ p = − µ ˜ p + β κ ˜ p → (cid:20) − ˜ p − ∂ ˜ p∂κ (cid:21) (24)˙˜ p → = [ βκ − µ ] ˜ p → − β κ ˜ p → (cid:18) ˜ p → + ∂ ˜ p → ∂κ (cid:19) (25)In order to solve these equations we need some ansatz for the dependence of the mean values ˜ p and ˜ p → on the meandegree in the graph. Inspired by the whole derivation of the Cavity Master Equation for spin systems, we propose˜ p ( t ) = 12 [1 + tanh ( β κ χ ( t ) − µ )] (26)˜ p → ( t ) = 12 [1 + tanh ( β κ ǫ ( t ) − µ )] (27)where χ ( t ) and ǫ ( t ) are some time-dependent fields that are obtained by inverting these very formulas in terms of ˜ p and ˜ p → .From this ansatz we can express the derivatives with respect to the degree in (24) and (25) in terms of ˜ p and ˜ p → ,respectively (see appendix B). We get then a closed system of differential equations for these probabilities, that canbe solved numerically. The results of the integration are shown in figure 6. VI. EXTENSIONS TO SIR-SIRS MODELS
Models with more compartments are also ubiquitous in epidemics modeling. In particular the SIR family, includingSEIR and SIRS, that have experienced a boost in attention due to the COVID-19 currently ongoing epidemic. TheCavity Master Equation, that was introduced for spins σ i = ± P i ( σ i )d t = X σ ′ X σ ∂i h ( − δσ i ,σ ′ r i ( σ ′ , σ ∂i ) (cid:2) Y k ∈ ∂i P ki ( σ k | σ ′ ) (cid:3) P i ( σ ′ ) i (28)d p ij ( σ i | σ j )d t = X σ ′ X σ ∂i \ j " ( − δσ i ,σ ′ r i [ σ ′ , σ ∂i ] h Y k ∈ ∂i \ j p ki ( σ k | σ ′ ) i p ij ( σ ′ | σ j ) (29)Sums run over all q states of each variable ( q ∈ { S, I, R } for SIR and SIRS), and the ( −
1) factor chooses betweenterms that contribute to the state σ ′ = σ i , and terms that take σ i to some other configuration σ ′ . For instance, inthe SIR and SIRS cases the equations reduce to (after considering the explicit form of the rate functions):d P i ( σ i ≡ I )d t = − µP i ( I ) + βP i ( S ) X k ∈ ∂ i p ki ( I | S ) (30)d P i ( S )d t = γP i ( R ) − βP i ( S ) X k ∈ ∂ i p ki ( I | S ) (31)d P i ( R )d t = − γP i ( R ) + µP i ( I ) (32)d p ij ( I | S )d t = − µp ij ( I | S ) + βp ij ( S | S ) X k ∈ ∂ i \ j p ki ( I | S ) (33)d p ij ( S | S )d t = γp ij ( R | S ) − βp ij ( S | S ) X k ∈ ∂ i \ j p ki ( I | S ) (34)d p ij ( R | S )d t = − γp ij ( R | S ) + µp ij ( I | S ) (35)Constants µ and β are still the recovery rate and the infection rate, while γ is the rate at which immunity is lostand patients pass from the recovered compartment back into the susceptible one. In the case of SIRS model γ > γ = 0.The first three equations (30,31,32) are still the same obtained by the dynamical message passing in [5] but, as inthe SIS case, the conditional distribution equation (33) differs from that in rDMPd p ij ( I | S )d t = − µp ij ( I | S ) + βP j ( S ) X k ∈ ∂ i \ j p ki ( I | S ) (36)in the factor P j ( S ) before the sum over neighbors. Notice that this difference implies that only four equations areneeded in the rDMP case: 3 for each node ˙ P i ( I ) , ˙ P i ( S ) , ˙ P i ( R ) and one on each directed edge ˙ p ij ( I | S ); while in ourcase all 6 equations need to be integrated (3 on nodes, 3 on edges) since they are all mutually dependent.Numerical experiments show that the CME approach fits better the average evolution of Monte Carlo simulations(figure 7) for SIRS and SIR models. VII. CONCLUSIONS
The analytical and numerical results on susceptible-infectious-susceptible model have shown that cavity masterequation (CME) is an effective approximation to the average dynamics of epidemic systems, where average is intendedover many stochastic realizations of the epidemics process. It outperforms individual and pair based mean fieldapproximations, and corrects a term in the dynamic message passing approach. We have explored with some successthe average case (topology) predictions for random regular and Erdos-Renyi graphs, and also showed that the approachis extensible to more elaborated models as SIR and SIRS.Many aspects remain to be studied. First, although not shown, we have seen that CME (and probably also allothers methods) fails when collective behavior appears, for instance, in the edge of the endemic transition. This isexpected since a key step of the cavity master equation is the assumption of some sort of factorization. What is more0 P r ob -I n f- tMC simCMErDMPIBMF 0 0.05 0.1 0.15 0.2 0 10 20 30 40 50 P r ob - R ec - t 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0 50 100 150 200 I n f ec t e d - s i ze t MC simCMEIBMF 0 0.2 0.4 0.6 0.8 1 0 50 100 150 200 R ec ov e r e d - s i ze t Fig. 7. Left : Probability of node 29 to be infected as function of time, discounting at each time the situations inwhich epidemic disappear in a SIRS model.
Right : Percentage of the nodes in the infected state for a SIR model.Both figures are for Zacharia’s graph, with β = 0 . γ = 0 .
05 and both inserted graphics show the correspondingpercentage of nodes in the recovered state.surprising is that the approach is particulary bad in 1D ring networks. Exploring why and how to improve in suchcases seems an appealing direction to us.Second, the average case solution is far from general, relying in some physical intuition to close the set of equations.The general case, in the graphs sense and also in the model sense (SIR, SIRS, SEIR), remains to be studied. Whetherwe can turn CME equations into some sort of jacobian stability approach for a general type of graphs seems also agood approach.Third, the extension of CME to Potts models with q states need to be formalized in the framework of point processesthat was used to derive CME in the first place. While doing this, attention should be put on the possibility of writinga CME-likelihood that could be used for inverse problems or for inference with partial information over some nodes,which is a key issue in current projects for testing and tracing.
Appendix A: Cavity master equation for SIS model.
In this appendix we will present more details on the formulation of the CME for the SIS model. Let’s take asstarting points equations (9) and (10):d P ( σ i )d t = − X σ ∂i h r i ( σ i , σ ∂i ) (cid:2) Y k ∈ ∂i P ( σ k | σ i ) (cid:3) P ( σ i ) − r i ( − σ i , σ ∂i ) (cid:2) Y k ∈ ∂i P ( σ k | − σ i ) (cid:3) P ( − σ i ) i (A1)d p ( σ i | σ j )d t = − X σ ∂i \ j " r i [ σ i , σ ∂i ] h Y k ∈ ∂i \ j p ( σ k | σ i ) i p ( σ i | σ j ) − r i [ − σ i , σ ∂i ] h Y k ∈ ∂i \ j p ( σ k | − σ i ) i p ( − σ i | σ j ) (A2)Due to the complementarity of the probabilities P i ( I ) and P i ( S ) we can focus on just deriving the equation for oneof the terms. Let’s take P i ( σ ≡ I ) and rewrite equation (9) accordingly.d P i ( I )d t = − X σ ∂i h r i ( I, σ ∂i ) (cid:2) Y k ∈ ∂i P ki ( σ | I ) (cid:3) P i ( I ) − r i ( S, σ ∂i ) (cid:2) Y k ∈ ∂i P ki ( σ k | S ) (cid:3) P i ( S ) i (A3)The rate r i ( I, σ ∂i ) is the transition rate from state I → S and in the SIS model and in most of epidemic models itdoes not depends on the state of the contacts of the Infected person. It just depends on the typical time of recovering(or die) from the disease, so this rate is directly the recovering rate r i ( I, σ ∂i ) = µ . On the other hand r i ( S, σ ∂i ) isthe transition rate from S → I . In infectious diseases, the contagion can only occur if a susceptible individual is incontact with an infected one ( r i ( S, σ k ≡ I ) = β , r i ( S, σ k ≡ S ) = 0), and the probability of getting infected is additivewith the number of infected contacts. This means that we can rewrite r i ( S, σ ∂i ) = P k ∈ ∂i r i ( S, σ k ).d P i ( I )d t = − µP i ( I ) X σ ∂i Y k ∈ ∂i P ki ( σ | I ) + P i ( S ) X σ ∂i h X k ′ ∈ ∂i r i ( S, σ k ′ ) i Y k ∈ ∂i P ki ( σ k | S ) (A4)1Let’s apply the equivalence P σ ∂i Q k ∈ ∂i → Q k ∈ ∂i . P σ k to invert the sum and product in the first term of the righthand side (RHS).d P i ( I )d t = − µP i ( I ) Y k ∈ ∂i h X σ k P ki ( σ | I ) i + P i ( S ) X σ ∂i h X k ′ ∈ ∂i r i ( S, σ k ′ ) i Y k ∈ ∂i P ki ( σ k | S ) (A5)The term P σ k P ki ( σ | I ) is exactly the sum of P ki ( I | I ) + P ki ( S | I ) = 1 and therefore also the product is equal to 1.d P i ( I )d t = − µP i ( I ) + P i ( S ) X σ ∂i h X k ′ ∈ ∂i r i ( S, σ k ′ ) i Y k ∈ ∂i P ki ( σ k | S ) (A6)Now we rewrite the second term of the RHS as follows:d P i ( I )d t = − µP i ( I ) + P i ( S ) X σ ∂i X k ′ ∈ ∂i r i ( S, σ k ′ ) P k ′ i ( σ | S ) Y k ∈ ∂i \ k ′ P ki ( σ | S ) (A7)d P i ( I )d t = − µP i ( I ) + P i ( S ) X σ k ′ X k ′ ∈ ∂i r i ( S, σ k ′ ) P k ′ i ( σ | S ) X σ ∂i \ k ′ Y k ∈ ∂i \ k ′ P ki ( σ | S ) (A8)Exchanging again the sum and the product, but now on the last term, we get P σ ∂i \ k ′ h Q k ∈ ∂i \ k ′ P ki ( σ | S ) (cid:3)i = 1, andusing that r i ( S, S ) = 0 and r i ( S, I ) = β we obtain equation (11):d p i d t = − µp i + β (1 − p i ) X k p ki (A9)Following the same steps is easy to derive the equation (12) for conditional probabilities. Appendix B: Average case equations for Erdos-Renyi graphs
In this appendix we will show how to perform the closure of average-case equations for Erdos-Renyi graphs (sub-section V B). Let’s start by explicitly writing the derivatives with respect to κ : ∂ ˜ p∂κ = ∂∂κ X γ P ( γ ) p γ ! = ∂∂κ X γ e − κ κ γ +1 ( γ + 1)! p γ ! (B1) ∂ ˜ p → ∂κ = ∂∂κ X γ g link ( γ ) p γ → ! = ∂∂κ X γ e − κ κ γ γ ! p γ → ! (B2)Computation of (B1) and (B2) gives: ∂ ˜ p∂κ = − X γ e − κ κ γ +1 ( γ + 1)! p γ + X γ ( γ + 1) e − κ κ γ ( γ + 1)! p γ = − ˜ p + 1 κ X γ ( γ + 1) P ( γ ) p γ (B3) ∂ ˜ p → ∂κ = − X γ e − κ κ γ ( γ )! p γ → + X γ γ e − κ κ γ − γ ! p γ → = − ˜ p → + 1 κ X γ γ g link ( γ ) p γ → (B4)The sums in the right hand sides of (B3) and (B4) are also involved in equations (20) and (21). Then, rememberingthat κ = P γ ( γ + 1) P ( γ ), we can rewrite equation (20) as follows:˙˜ p = − µ ˜ p + β ˜ p → X γ ( γ + 1) P ( γ )(1 − p γ )˙˜ p = − µ ˜ p + β ˜ p → X γ ( γ + 1) P ( γ ) − β ˜ p → X γ ( γ + 1) P ( γ ) p γ ˙˜ p = − µ ˜ p + βκ ˜ p → − βκ ˜ p → (cid:18) ˜ p + ∂ ˜ p∂κ (cid:19) (B5)2which leads directly to equation (24). Equation (25) can be derived by an analogous procedure, using an equivalentexpression for the mean connectivity: κ = P γ γ g link ( γ )Now we just need closed expressions for the derivatives in (24) and (25). In order to do so, let’s compute the κ -derivative on both sides of ansatz (26) and (27). We get: ∂ ˜ p∂κ = 12 (cid:2) − tanh ( β κ χ ( t ) − µ ) (cid:3) β χ ( t ) (B6) ∂ ˜ p → ∂k = 12 (cid:2) − tanh ( β κ ǫ ( t ) − µ ) (cid:3) βǫ ( t ) (B7)We can re-use equations (26) and (27) for eliminating χ ( t ) and ǫ ( t ) from (B6) and (B7), thus obtaining the followingclosed expressions for the derivatives: ∂ ˜ p∂κ = 12 κ h − (2˜ p − i (cid:2) tanh − (2˜ p −
1) + µ (cid:3) (B8) ∂ ˜ p → ∂k = 12 κ h − (2˜ p → − i (cid:2) tanh − (2˜ p → −
1) + µ (cid:3) (B9)This allows to numerically solve equations (24) and (25). [1] William Ogilvy Kermack and Anderson G McKendrick. A contribution to the mathematical theory of epidemics. Proceed-ings of the royal society of london. Series A, Containing papers of a mathematical and physical character , 115(772):700–721,1927.[2] Péter L Simon, Michael Taylor, and Istvan Z Kiss. Exact epidemic models on graphs using graph-automorphism drivenlumping.
Journal of mathematical biology , 62(4):479–508, 2011.[3] Faryad Darabi Sahneh, Caterina Scoglio, and Piet Van Mieghem. Generalized epidemic mean-field model for spreadingprocesses over multilayer complex networks.
IEEE/ACM Transactions on Networking , 21(5):1609–1620, 2013.[4] Andrey Y Lokhov, Marc Mézard, Hiroki Ohta, and Lenka Zdeborová. Inferring the origin of an epidemic with a dynamicmessage-passing algorithm.
Physical Review E , 90(1):012801, 2014.[5] Munik Shrestha, Samuel V Scarpino, and Cristopher Moore. Message-passing approach for recurrent-state epidemic modelson networks.
Physical Review E , 92(2):022821, 2015.[6] Fabrizio Altarelli, Alfredo Braunstein, Luca Dall’Asta, Alejandro Lage-Castellanos, and Riccardo Zecchina. Bayesianinference of epidemics on networks via belief propagation.
Physical review letters , 112(11):118701, 2014.[7] Alfredo Braunstein and Alessandro Ingrosso. Inference of causality in epidemics on temporal contact networks.
Scientificreports , 6:27538, 2016.[8] Werner Krauth and Marc Mézard. The cavity method and the travelling-salesman problem.
EPL (Europhysics Letters) ,8(3):213, 1989.[9] Erik Aurell, Gino Del Ferraro, Eduardo Domínguez, and Roberto Mulet. Cavity master equation for the continuous timedynamics of discrete-spin models.
Physical Review E , 95(5):052119, 2017.[10] Eric Cator and Piet Van Mieghem. Second-order mean-field susceptible-infected-susceptible epidemic threshold.
Physicalreview E , 85(5):056111, 2012.[11] Kieran J Sharkey. Deterministic epidemiological models at the individual level.