Memory-induced complex contagion in epidemic spreading

The miccSIS model describes a population of agents that can be either susceptible (healthy) or infected (infectious), and is embedded on an undirected, unweighted network. The contact network is encoded in the adjacency matrix, which is nonspatial (it carries no information about the agents' physical position) and static (it remains fixed over time).

Infected nodes have a constant infectivity rate, υ, and continuously spread doses of contagion towards their entire neighborhood. They target all of their neighbors equally, transmitting pathogen along each edge at constant rate υ. Susceptible nodes collect these toxins from all their neighbors, amassing a total viral load κ, and transition to the infected state with probability ${\psi }_{\inf }^{* }(\kappa ){\rm{d}}\kappa$, where ${\psi }_{\inf }^{* }(\kappa )$ is the infection probability density. Infected nodes are unaffected by the toxins (their viral load does not increase) and recover spontaneously after a random time t, with interevent time distribution ${\psi }_{\mathrm{rec}}(t)$. At recovery, their viral load is completely erased and they become susceptible once again.

As in the standard SIS model, susceptible nodes whose nearest neighborhood is completely healthy cannot become infected. Since no active processes are associated to their state, they are irrelevant for the immediate evolution of the system. However, these inactive nodes play a crucial role in the long-term dynamics of the miccSIS model, therefore we assign them to an additional compartment, which we call dormant. A dormant node transitions to susceptible as soon as one of its neighbors becomes infected. Conversely, when the last infected neighbor of a susceptible node recovers, the latter transitions to the dormant state. At this point, the viral load it had previously amassed starts to deteriorate with relaxation time ζ, modeling its long-term memory. This feature mimics the restoring of an individual's immune system, or the gradual loss of interest of an opinion, idea, or trend.

In summary, infected (I) agents spread pathogen to all their neighbors and recover spontaneously. While susceptible (S) agents have at least one infected neighbor and continuously accumulate viral load, dormant (D) agents have a fully healthy neighborhood and cannot become infected. There are two types of active processes which entail one or possibly more transitions (see figure 1):

• Infection of susceptible agent j. Agent j transitions from susceptible to infected. Additionally, all of j's neighbors that were dormant transition to susceptible (and resume their accumulation of viral load).
• Recovery of infected agent j. If all of j's neighbors are healthy, j transitions from infected to dormant. If at least one of j's neighbors is infected, j transitions from infected to susceptible. Additionally, all of j's neighbors that were susceptible and had only one infected neighbor (i.e. agent j) transition to dormant (and their viral load starts to decay).

Infected agents are unaffected by the viral load, and ignore any new doses received from their infected neighbors. When an infected agent recovers it erases all the previously amassed viral load.

Zoom In Zoom Out Reset image size

Overall, the system's evolution is determined by a set of discrete, stochastic processes, an infection for each susceptible node and a recovery for each infected node. All these processes are statistically independent, which enables the use of the generalized non-Markovian Gillespie algorithm [27], capable of simulating memoryful dynamics in continuous time. A key ingredient of this algorithm is the instantaneous hazard rate of an active process, ω(t), obtained from its interevent time distribution, ψ(t), and corresponding survival probability, ${\rm{\Psi }}(t)={\int }_{t}^{\infty }\psi (t^{\prime} ){\rm{d}}t^{\prime}$, as ω(t) = ψ(t)/Ψ(t). In short, ω(t) measures the probability per unit of time that the corresponding event takes place between t and $t+{\rm{d}}t$ [36]. For a Poisson point process, the interevent time distribution is exponential and the hazard rate is, therefore, constant. In general, interevent time distributions decaying slower (respectively, faster) than exponential lead to asymptotically decreasing (increasing) hazard rates.

While recoveries are readily incorporated into this framework, ${\omega }_{\mathrm{rec}}(t)={\psi }_{\mathrm{rec}}(t)/{{\rm{\Psi }}}_{\mathrm{rec}}(t)$, infections require some additional attention. Since the activity in a susceptible node's neighborhood varies over time, its instantaneous amassment rate, $\tilde{\upsilon }(t)$, is not constant. Then, the instantaneous hazard rate for infections is

$$\begin{eqnarray}&&{\omega }_{\inf }(t)=\tilde{\upsilon }(t)\displaystyle \frac{{\psi }_{\inf }^{* }(\kappa (t))}{{{\rm{\Psi }}}_{\inf }^{* }(\kappa (t))},\end{eqnarray} \tag{ 1 }$$

with κ(t) the accumulated viral load at time t (see appendix A for a detailed derivation).

In general, the infectivity rate, υ, the relaxation time, ζ, the infection probability density, ${\psi }_{\inf }^{* }$, and the recovery interevent time distribution, ψ_rec, may vary from node to node. For example, one could model distinct age groups by segregating the population and assigning different values of the parameters to each subpopulation. Notwithstanding, in order to eliminate the effects of node heterogeneities, in the present work we use the same υ, ζ, ${\psi }_{\inf }^{* }$, and ${\psi }_{\mathrm{rec}}$ for all nodes. For this same reason, we limit our analysis to random degree-regular networks, particularly with degree k = 4.

For infections we select the versatile Weibull distribution, with shape parameter α and scale parameter μ

$$\begin{eqnarray}&&{\psi }_{\inf }^{* }(\kappa )=\alpha {\mu }^{\alpha }{\kappa }^{\alpha -1}{{\rm{e}}}^{-{\left(\mu \kappa \right)}^{\alpha }}.\end{eqnarray} \tag{ 2 }$$

When α > 1 it presents a peak, resembling a bell curve, α = 1 corresponds to a Poisson distribution, and for α < 1 it has power-law-like fat tails. The corresponding instantaneous hazard rate is

$$\begin{eqnarray}&&{\omega }_{\inf }(t)=\upsilon \alpha {\mu }^{\alpha }z(t){\left[\kappa (t)\right]}^{\alpha -1},\end{eqnarray} \tag{ 3 }$$

with z(t) the number of infected neighbors at time t (see appendix A for a detailed derivation). With α > 1 (respectively, α < 1), ω_inf(t) increases (decreases) monotonically with κ(t).

Notice that when α = 1 we recover the customary expression of the standard SIS model, ${\omega }_{\inf }(t)\sim z(t)$. For $\alpha \ne 1$, however, equation (3) cannot be written as a linear superposition of independent transmission channels. In the particular case of α = 2, for example, the hazard rate includes a quadratic term, ${\omega }_{\inf }(t)={az}(t)+b{\left[z(t)\right]}^{2}$. Hence, the agents' memory induces a complex contagion scheme, even though the model does not explicitly incorporate any social reinforcement/inhibition mechanisms (see appendix B for a detailed discussion). To further isolate the effects of the modified infection mechanism, we treat recoveries as Poisson processes with rate η, i.e. with constant hazard rate ω_rec(t) = η.

We define the effective spreading ratio, λ, as the average time required to recover over the average viral load needed to become infected, nondimensionalized by the infectivity rate, $\lambda =\upsilon \langle {t}_{\mathrm{rec}}\rangle /\langle {\kappa }_{\inf }\rangle$. With the selected distributions we find an expression for the scale parameter, $\mu =\eta \lambda {\left(\alpha \upsilon \right)}^{-1}{\rm{\Gamma }}({\alpha }^{-1})$, with Γ the gamma function. Notice that, once again, α = 1 recovers the customary expression of the standard SIS model, $\lambda =\upsilon \mu /\eta$ (with infectious rate $\upsilon \mu$ per infected neighbor). As for the relaxation time of dormant nodes' viral load, we consider the limit cases of instantaneous decay, ζ = 0, and perpetual accumulation, $\zeta =\infty$.

For illustrative purposes, consider the system depicted in figure 2(a), where all nodes are initially healthy except for D. Node C becomes infected at time t₀ and subsequently infects B at t₁. During the interval $t\in [{t}_{0},{t}_{1}]$, node A's viral load, κ_A, grows with rate υ, but from t₁ onwards it will increase with rate 2υ. At t₂, node C recovers and κ_A reduces its accumulation rate back to υ, and when B recovers at t₃, κ_A starts to decay with relaxation time ζ. Finally, C becomes infected once again at t₄ and κ_A resumes its growth at rate υ. Figures 2(b)–(e) show the evolution of κ_A and ω_A for various ζ and α. Hereafter we use temporal units such that η = 1 and, without loss of generality, set $\upsilon =1$.

Zoom In Zoom Out Reset image size

We begin our analysis for ζ = 0, with dormant nodes instantly erasing their viral load. Thus, when the outbreak reenters their neighborhood, they become susceptible starting afresh (with κ = 0). Hence, the only memory effect present is during the infection period, which we interpret as a short-term memory mode.

Each node is uniquely defined by its state, infected or healthy, and its accumulated viral load. The state of a node only changes with the transitions (i) infected to healthy, and (ii) susceptible to infected. When a node transitions between susceptible and dormant, its state remains unaltered. On the other hand, the viral load (i) is erased instantly when an infected node recovers, (ii) increases proportionally to the number of infected neighbors while a node is susceptible, and (iii) is erased instantly when a susceptible node becomes dormant.

In the supplemental material (SM), we derive exact stochastic evolution equations for both the state and accumulated viral load of each node. Applying a mean-field approximation, these equations can be reduced to a single expression for the late-time prevalence, ρ, the average fraction of nodes that are infected in the steady-state (see SM for a detailed derivation). If all nodes have the same degree k, this equation reads

$$\begin{eqnarray}&&-\rho +{ak}\rho (1-\rho )-{bk}\rho {\left(1-\rho \right)}^{k}=0,\end{eqnarray} \tag{ 4 }$$

with a > 0. Roughly speaking, the first term corresponds to the recovery of infected nodes, the second term to the infection of susceptible nodes, and the third is related to susceptible nodes becoming dormant. The shape factor of the infection probability, α, and the effective spreading ratio, λ, are encapsulated in a and b, which become constant parameters when ρ ≈ 0.

Equation (4) uncovers the existence of an epidemic threshold, where the population transitions from a healthy, pathogen-free state to an endemic, disease burdened one. Linear stability analysis reveals that the healthy phase, ρ = 0, becomes unstable at $a-b-1/k=0$, yielding a phase transition at the critical point ${a}_{{\rm{c}}}=b+1/k$. Moreover, the nature of the transition changes at $a-{bk}=0$, yielding a tricritical point at ${a}_{\mathrm{tc}}=1/(k-1)$, ${b}_{\mathrm{tc}}=1/k(k-1)$. Then the steady state is endemic for a > a_c, and the phase transition is continuous for b < b_tc and discontinuous for b > b_tc. Finally, the prevalence of the endemic phase scales as $\rho \propto {\left(a-{a}_{{\rm{c}}}\right)}^{\beta }$, with β_c = 1 when the transition is continuous, and β_tc = 1/2 at the tricritical point.

Figure 3 shows the phase diagram in terms of the original parameters α and λ. We observe that the critical point initially increases monotonically with α, but afterwards saturates for $\alpha \to \infty$. This result is consistent with the possibly largest epidemic threshold reported in [37]. On the other hand, the transition to endemicity is discontinuous for α > α_tc ≈ 2.348. Recall that the infection probability density presents a peak for α > 1, and, in fact, tends towards a delta function in the limit $\alpha \to \infty$. Thus, a node requires a quasi deterministic amount of viral load to become infected, mimicking a threshold model, which commonly exhibits a discontinuous phase transition [38].

Zoom In Zoom Out Reset image size

In order to verify these analytic findings, we perform extensive stochastic simulations. We begin by exploring the position of the critical point, ${\lambda }_{{\rm{c}}}$, that separates an absorbing (healthy) phase (λ < λ_c) from an active (endemic) one (λ > λ_c). The simulations start well into the active phase with a fully infected population, and quasistatically decrease the control parameter, λ, until finite-size fluctuations trap the system in the absorbing state. We sample the late-time prevalence, ${\rho }_{\mathrm{st}}={\mathrm{lim}}_{t\to \infty }{N}_{{\rm{I}}}(t)/N$, of 10⁴ states, time-averaged over various trajectories (see SM for simulation details).

As shown in figure 4, λ_c indeed increases with α, and saturates for large values of α. Moreover, for α < 1 the approach to the critical point is very similar to the standard SIS (α = 1), consistent with a continuous phase transition with exponent β = 1. On the other hand, for α = 4 the curves terminate at a remarkably high prevalence, consistent with a discontinuous phase transition. Finally, the curves for α = 2 also terminate at a rather high prevalence, which deviates from the analytic prediction. Nonetheless, since the apparent discontinuity decreases with the system size, this observation is most likely related to finite-size effects.

Zoom In Zoom Out Reset image size

We complement these results with the analysis of patient zero scenarios, the arrival of an infected agent in a previously unaffected population. We employ the lifespan method [39], which simulates outbreaks starting from a single infected node. Each single-seed realization is characterized by its coverage, K, defined as the number of distinct nodes that have become infected at least once, and can be either finite or endemic. While the former return to the absorbing state, the latter evolve towards an active steady state. In finite systems, we introduce a coverage threshold, ${K}_{\mathrm{th}}={c}_{\mathrm{th}}N$, with 0 < c_th < 1. A realization is declared endemic whenever its coverage reaches the threshold, and those that terminate without surpassing it are considered finite. As reported in [39], the value of c_th does not modify the qualitative results. Hereafter, we use c_th = 0.75.

For a fixed value of λ we run 10⁴ realizations, each starting with a single, randomly chosen infected node, and a system cleared of all viral load. If an outbreak becomes endemic, we extend the simulation until it reaches the steady-state, and then measure the late-time prevalence, ${\rho }_{\mathrm{st}}^{* }$ (see SM for simulation details). The results are shown in figures 5(a), (b), which include the previously computed ρ_st. While ${\rho }_{\mathrm{st}}^{* }$ shows a discontinuous jump for α = 4, with α = 2 it grows continuously and coincides with ρ_st, supporting our analytical findings. Moreover, since α = 2 is close to tricritical point, a cross-over effect towards the exponent β_tc is expected, explaining the rather steep approach towards the critical point.

Zoom In Zoom Out Reset image size

Overall, these analytic and simulated results indicate that the system's macroscopic properties are drastically affected by the microscopic details of the infection mechanism. In particular, the critical point that separates the healthy from the endemic phase grows with α, and the nature of the phase transition changes from continuous to discontinuous at the tricrital point α_tc.

We first analyze α > 1, for which the infection probability presents a peak and the instantaneous infection rate increases monotonically with the accumulated viral load. The patient zero results are shown in figures 7(a), (b), which include the previously computed ρ_st, and the late-time prevalence of single-seed outbreaks that are able to become endemic, ρ_st* (see SM for simulation details). We find that the average coverage, $\bar{c}$, and the endemic probability, P₁, coincide for all values of λ and present a continuous phase transition at ${\lambda }_{{\rm{c}}}(\bar{c})$. However, this point is notably larger than the critical point of the late-time prevalence, λ_c(ρ_st). While ρ_st presents a continuous phase transition, ${\rho }_{\mathrm{st}}^{* }$ exhibits a discontinuous phase transition at ${\lambda }_{{\rm{c}}}({\rho }_{\mathrm{st}}^{* })={\lambda }_{{\rm{c}}}(\bar{c})$. As expected, the two prevalence curves overlap after the abrupt jump.

Zoom In Zoom Out Reset image size

This evidences the existence of an intermediate region $\lambda \in [{\lambda }_{{\rm{c}}}({\rho }_{\mathrm{st}}),{\lambda }_{{\rm{c}}}(\bar{c})]$ where all single-seed outbreaks return to the absorbing state, whereas fully infected populations evolve towards an active steady state. The key ingredient to understand this phenomenon is the environment of frozen viral load. During the simulations that measure the late-time prevalence, the viral loads are well thermalized, enabling the outbreak to remain in an active state. Conversely, this environment is deficient in single-seed outbreaks, as the system has not yet reached its steady state. Hence, outbreaks are unable to produce sufficient new infections and rapidly become trapped in the absorbing state.

These results indicate that the system displays two attractors in this intermediate region. Then, for $\zeta =\infty$ and α > 1 the system's phase diagram exhibits an additional bistable phase, that separates the usual healthy and endemic phases. The associated hysteresis loop, however, has a rather exotic nature: although its lower branch presents the expected discontinuity, the upper branch connects the two attractors in a continuous manner. This contrasts with previous findings of bistability, were the hysteresis loop is bounded by two discontinuities [24, 29, 40]. Moreover, the transition to full endemicity is hybrid [41]: the endemic probability grows continuously at ${\lambda }_{{\rm{c}}}({\rho }_{\mathrm{st}}^{* })$, but the late-time prevalence jumps discontinuously.

Finally, we study α < 1, for which the infection probability presents power-law-like fat tails an the instantaneous infection rate decreases monotonically with the accumulated viral load. In figures 8(a), (b) we show the patient zero analysis for the standard SIS (α = 1) and broad-tailed infection distributions (α < 1). Here, we additionally compute P₃, the probability that a single-seed outbreak reaches the coverage threshold three times (see SM for simulation details). With α = 1, P₁ and P₃ are practically identical, indicating that an outbreak that surpasses the coverage threshold once remains active long enough to surpass the threshold two more times. Thus P₁ is an adequate proxy for the true endemic probability, P_∞, which additionally coincides with $\bar{c}$ and ρ_st.

Zoom In Zoom Out Reset image size

For α < 1 the situation is quite different. Firstly, the average coverage starts growing continuously at ${\lambda }_{{\rm{c}}}(\bar{c})$, when all other order parameters are still identically zero. Additionally, the transition point of P₁ is significantly lower than that of P₃. Thus, there is a wide interval where all outbreaks that surpass the threshold once eventually terminate in the absorbing state, evidencing the inadequateness of P₁ as a measure of the true endemic probability. The inflection point of P₃ is much closer to the transition point of ρ_st, which suggests that the critical point of the endemic probability (${P}_{\infty }$, the probability to surpass the threshold an infinite amount of times) coincides with ${\lambda }_{{\rm{c}}}({\rho }_{\mathrm{st}})$. Beyond this point, P_∞ is expected to coincide with $\bar{c}$, indicating that the endemic probability presents a discontinuous phase transition. Then, the transition to full endemicity is again hybrid. Nonetheless, here the late-time prevalence grows continuously, while the endemic probability jumps discontinuously.

In this case, we find an intermediate region $\lambda \in [{\lambda }_{{\rm{c}}}(\bar{c}),{\lambda }_{{\rm{c}}}({\rho }_{\mathrm{st}})]$ where outbreaks are unable to become endemic (${P}_{\infty }=0$) but affect a macroscopic fraction of the population (c > 0). In figures 8(c), (d) we show the prevalence, $\bar{\rho }(t)=\langle {N}_{{\rm{I}}}(t)\rangle /N$, of realizations that surpass the coverage threshold once, for λ = 0.33 and λ = 0.4. We include the results for ζ = 0 for comparison (see SM for simulation details). Both values of λ are located in the active phase of the short-term memory mode (ζ = 0), and so the endemic realizations converge monotonically towards their active steady states. For the long-term memory mode ($\zeta =\infty$), λ = 0.33 is located in the intermediate region. In this case, we observe that outbreaks grow up to a maximum, after which their prevalence gradually diminishes until they reach the absorbing state. This behavior is typically observed in SIR-like dynamics and is reminiscent of excitable media. In the endemic phase (λ = 0.4) the outbreaks continue presenting a peak, but afterwards relax towards an active steady state.

In conclusion, for $\zeta =\infty$ and α < 1 the usual healthy and endemic phases are separated by an additional excitable phase. This excitable behavior is again a consequence of the environment of frozen viral load. Independently of ζ, an outbreak starts from a single infected node in a population cleared of viral load. Then it initially evolves as if the agents only had the short-term memory mode (clearly appreciable for $t\in [0,25]$ in figure 8(c)), rapidly achieving a large coverage. When the outbreak revisits a previously affected area, the long-term memory mode is activated and the frozen viral load impedes new infections. Thus, dormant nodes are effectively removed from the dynamic, impede the outbreak to grow, and eventually cause its extinction. To the extend of our knowledge, this excitable behavior has not been previously reported in comparable SIS-like models.

Here, we summarize the derivation reported in [27]. Consider a set of M statistically independent, discrete, stochastic processes, each with an interevent time distribution ψ_j(t) and corresponding survival probability ${{\rm{\Psi }}}_{j}(t)={\int }_{t}^{\infty }{\psi }_{j}(t^{\prime} ){\rm{d}}t^{\prime}$. At a certain moment in time t₀, process j has been active for t_j units of time. Let $\phi (\tau ,i| \left\{{t}_{k}\right\}){\rm{d}}\tau$ denote the joint probability that the next-occurring event takes place in the interval $t\in [{t}_{0}+\tau ,{t}_{0}+\tau +{\rm{d}}\tau ]$ and corresponds to process i, conditioned by the set of elapsed times $\left\{{t}_{k}\right\}$. This probability density can be expressed as

$$\begin{eqnarray}&&\phi (\tau ,i| \left\{{t}_{k}\right\})=\displaystyle \frac{{\psi }_{i}({t}_{i}+\tau )}{{{\rm{\Psi }}}_{i}({t}_{i}+\tau )}{\rm{\Phi }}(\tau | \left\{{t}_{k}\right\}),\end{eqnarray} \tag{ A1 }$$

where

$$\begin{eqnarray}&&{\rm{\Phi }}(\tau | \left\{{t}_{k}\right\})=\displaystyle \prod _{j=1}^{M}\displaystyle \frac{{{\rm{\Psi }}}_{j}({t}_{j}+\tau )}{{{\rm{\Psi }}}_{j}({t}_{j})}\end{eqnarray} \tag{ A2 }$$

is the survival probability of τ, i.e. the conditional probability that no event takes place before t₀ + τ. Then the probability that the next event takes place in the interval $t\in [{t}_{0},{t}_{0}+\tau ]$ is

$$\begin{eqnarray}&&{\rm{\Xi }}(\tau | \left\{{t}_{k}\right\})=1-{\rm{\Phi }}(\tau | \left\{{t}_{k}\right\}).\end{eqnarray} \tag{ A3 }$$

Once the interval τ is known, the probability that the next-occurring event corresponds to process i is given by

$$\begin{eqnarray}&&{\rm{\Pi }}(i| \tau ,\left\{{t}_{k}\right\})=\displaystyle \frac{{\omega }_{i}({t}_{i}+\tau )}{{\displaystyle \sum }_{j=1}^{M}{\omega }_{j}({t}_{j}+\tau )},\end{eqnarray} \tag{ A4 }$$

with ${\omega }_{j}(t)={\psi }_{j}(t)/{{\rm{\Psi }}}_{j}(t)$ the instantaneous hazard rate of process j.

Equations (A3) and (A4) provide an algorithm that generates statistically correct sequences of events: (i) draw the interval by solving ${\rm{\Xi }}(\tau | \left\{{t}_{k}\right\})=u$, with $u\in U(0,1)$, (ii) increase the system time as $t\leftarrow t+\tau$, (iii) draw the process from the discrete distribution ${\rm{\Pi }}(i| \tau ,\left\{{t}_{k}\right\})$, (iv) revise the list of active processes, and (v) update the set of elapsed times as ${t}_{j}\leftarrow {t}_{j}+\tau$ (setting t_j = 0 for newly activated processes).

Recoveries are straightforwardly incorporated into this framework, with the elapsed time t_j measuring the period since agent j became infected (i.e. this occurred at t = t₀ − t_j). On the other hand, infection processes require the translation of infection probability densities into interevent time distributions. Consider susceptible agent j, characterized by its infection probability density ${\psi }_{j}^{* }(\kappa )$ and corresponding survival probability ${{\rm{\Psi }}}_{j}^{* }(\kappa )={\int }_{\kappa }^{\infty }{\psi }_{j}^{* }(\kappa ^{\prime} ){\rm{d}}\kappa ^{\prime}$. Its interevent time distribution, ${\psi }_{j}(t)$, is given by the normalization condition

$$\begin{eqnarray}&&{\psi }_{j}(t){\rm{d}}t={\psi }_{j}^{* }(\kappa ){\rm{d}}\kappa .\end{eqnarray} \tag{ A5 }$$

Since the activity in j's neighborhood may vary over time, the rate at which it amasses viral load is generally nonconstant. At time t it has amassed κ_j(t) units of viral load and has z_j(t) infected neighbors. If the system remains unaltered in an interval dt, node j will amass an additional ${\rm{d}}\kappa ={\tilde{\upsilon }}_{j}{\rm{d}}t$, with ${\tilde{\upsilon }}_{j}(t)={\sum }_{i=1}^{{z}_{j}(t)}{\upsilon }_{i}$ its instantaneous amassment rate. Substituting in equation (A5) we find

$$\begin{eqnarray}&&{\psi }_{j}(t)={\tilde{\upsilon }}_{j}(t){\psi }_{j}^{* }({\kappa }_{j}(t)).\end{eqnarray} \tag{ A6 }$$

For the survival probability we have

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{j}(t)={\int }_{t}^{\infty }{\psi }_{j}(t^{\prime} ){\rm{d}}t^{\prime} ={\int }_{{\kappa }_{j}(t)}^{\infty }{\psi }_{j}^{* }(\kappa ^{\prime} ){\rm{d}}\kappa ^{\prime} ={{\rm{\Psi }}}_{j}^{* }({\kappa }_{j}(t)),\end{eqnarray} \tag{ A7 }$$

which yields its instantaneous hazard rate

$$\begin{eqnarray}&&{\omega }_{j}(t)={\tilde{\upsilon }}_{j}(t)\displaystyle \frac{{\psi }_{j}^{* }({\kappa }_{j}(t))}{{{\rm{\Psi }}}_{j}^{* }({\kappa }_{j}(t))}.\end{eqnarray} \tag{ A8 }$$

Note that we can always write $t={t}_{0}+\tau$, with t₀ the time at which the system was last updated, and τ ≥ 0. Then, the instantaneous amassment rate remains constant in the interval $[{t}_{0},t]$, ${\tilde{\upsilon }}_{j}(t)={\tilde{\upsilon }}_{j}({t}_{0})$, and ${\kappa }_{j}(t)={\kappa }_{j}({t}_{0})+\tau {\tilde{\upsilon }}_{j}({t}_{0})$.

In our work, infections are governed by a Weibull distribution with shape parameter α and scale parameter μ

$$\begin{eqnarray}&&{\psi }_{\inf }^{* }(\kappa )=\alpha {\mu }^{\alpha }{\kappa }^{\alpha -1}{{\rm{e}}}^{-{\left(\mu \kappa \right)}^{\alpha }}\quad {{\rm{\Psi }}}_{\inf }^{* }(\kappa )={{\rm{e}}}^{-{\left(\mu \kappa \right)}^{\alpha }},\end{eqnarray} \tag{ A9 }$$

and recoveries by a Poisson process with rate η

$$\begin{eqnarray}&&{\psi }_{\mathrm{rec}}(t)=\eta {{\rm{e}}}^{-\eta t}\quad {{\rm{\Psi }}}_{\mathrm{rec}}(t)={{\rm{e}}}^{-\eta t}.\end{eqnarray} \tag{ A10 }$$

Since all infectors have the same infectivity rate υ, the instantaneous amassment rate of a susceptible node j becomes ${\tilde{\upsilon }}_{j}(t)=\upsilon {z}_{j}(t)$, with z_j(t) the number of its neighbors that are infected at time t. So the instantaneous hazard rates for infections and recoveries are, respectively, ${\omega }_{\inf }(t)=\upsilon \alpha {\mu }^{\alpha }z(t){\left[\kappa (t)\right]}^{\alpha -1}$ and ω_rec(t) = η.

Simple contagion describes purely dyadic interactions, thus we can identify each edge that connects a healthy node with an infected one as an isolated transmission channel. Consider at time t a susceptible node j and its infected neighbor i, which became infected at t_i < t. The probability that node i infects node j within the interval $(t,t+{\rm{d}}t)$ is ${\omega }_{i\to j}(t| {t}_{i}){\rm{d}}t$, regardless of the rest of the system. If node j has z_j(t) infectors at time t, the previous statement holds for each of them.

The total probability that node j becomes infected at time t depends on all of its incoming transmission channels. Since these are statistically independent, we can write

$$\begin{eqnarray}&&{\omega }_{j}(t)=\displaystyle \sum _{i=1}^{{z}_{j}(t)}{\omega }_{i\to j}(t| {t}_{i}),\end{eqnarray} \tag{ B1 }$$

where ${\omega }_{j}(t)$ is the instantaneous hazard rate of node j's infection process (i.e. the probability per unit of time that node j becomes infected at time t). Using

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{j}(t)=\displaystyle \frac{1}{{z}_{j}(t)}\displaystyle \sum _{i=1}^{{z}_{j}(t)}{\omega }_{i\to j}(t| {t}_{i})\end{eqnarray} \tag{ B2 }$$

we can write equation (B1) as

$$\begin{eqnarray}&&{\omega }_{j}(t)={{\rm{\Omega }}}_{j}(t){z}_{j}(t),\end{eqnarray} \tag{ B3 }$$

thus the total hazard rate is proportional to the number of current infectors. If ${\omega }_{i\to j}(t| {t}_{i})=\beta$ are constants (and homogenous for all pairs of nodes), we recover the standard SIS model with the familiar expression ${\omega }_{j}(t)=\beta {z}_{j}(t)$. When the transmission rates ${\omega }_{i\to j}$ are time-dependent, the dynamics has memory effects; thus, simple contagion can be non-Markovian (as in [23], for example).

When the dynamics are described by interactions that are not strictly dyadic, the contagion becomes complex. These processes usually incorporate an explicit social reinforcement or inhibition mechanism. Although the classification of complex contagion processes is yet to be formalized, they can be broadly categorized into two groups:

• In edge-centric approaches, one still considers the transmission channel from infected node i to susceptible node j. Now, however, the transmission rate ${\omega }_{i\to j}$ is affected by the neighborhoods of i and/or j (for specific examples, see [28, 29, 44]). Considering only nearest-neighbors, the transmission rate from node i to node j at time t, ${\omega }_{i\to j}(t| {z}_{i}(t),{z}_{j}(t))$, is a function of their current infected neighbors, z_i(t) and z_j(t). Although we can still write the total hazard rate ω_j(t) as in equation (B1), the instantaneous average defined in equation (B2) has an explicit dependence on z_i(t) and z_j(t). Consequently, ω_j(t) can be superlinear (reinforcement) or sublinear (inhibition) with the number of current infectors, z_j(t).
• On the other hand, node-centric approaches forgo the notion of transmission channels and directly prescribe the instantaneous hazard rate, ${\omega }_{j}(t)$. These usually incorporate thresholds, such as ${\omega }_{j}(t)=\delta ({T}_{j}-{z}_{j}(t))$ [30] or ${\omega }_{j}(t)={z}_{j}(t){\rm{\Theta }}({T}_{j}-{z}_{j}(t))+\beta {\rm{\Theta }}({z}_{j}(t)-{T}_{j})$ [45], which explicitly evidence the nonlinearity of ω_j(t) with z_j(t).

Consider an isolated pair of nodes i and j in the miccSIS model. Both are healthy when node i becomes infected at time t_i. The total amount of viral load that j has amassed at time t > t_i is ${\kappa }_{j}(t)=\upsilon (t-{t}_{i})$. The instantaneous hazard rate of node j's infection is

$$\begin{eqnarray}&&{\omega }_{j}(t)=\upsilon \alpha {\mu }^{\alpha }{z}_{j}(t){\left[{\kappa }_{j}(t)\right]}^{\alpha -1},\end{eqnarray} \tag{ B4 }$$

but since j has only one infected neighbor, it is given solely by the exposure to node i: ${\omega }_{j}(t)={\omega }_{i\to j}(t| {t}_{i})$, with

$$\begin{eqnarray}&&{\omega }_{i\to j}(t| {t}_{i})=\upsilon \alpha {\mu }^{\alpha }{\left[\upsilon (t-{t}_{i})\right]}^{\alpha -1}.\end{eqnarray} \tag{ B5 }$$

Now consider an infected node j that at time t₀ recovers and becomes dormant (all its neighbors are healthy). At time t it has a set of $\left\{i\right\}$ infected neighbors that became infected at times $\left\{{t}_{i}\right\}$, with ${t}_{0}\lt {t}_{i}\lt t$, $\forall i=1,2,\,\ldots ,\,{z}_{j}(t)$. The total amount of viral load that j has amassed at time t is

$$\begin{eqnarray}&&{\kappa }_{j}(t)=\chi (t)+\displaystyle \sum _{i=1}^{{z}_{j}(t)}\upsilon (t-{t}_{i}),\end{eqnarray} \tag{ B6 }$$

where χ(t) stores the viral load accumulated from neighbors that became infected in the interval (t₀, t) but are currently healthy. Substituing equations (B6) in (B4), we find

$$\begin{eqnarray}&&{\omega }_{j}(t)=\upsilon \alpha {\mu }^{\alpha }{z}_{j}(t){\left[\chi (t)+\displaystyle \sum _{i=1}^{{z}_{j}(t)}\upsilon (t-{t}_{i})\right]}^{\alpha -1},\end{eqnarray} \tag{ B7 }$$

which, using equation (B5), can be written for $\alpha \ne 1$ as

$$\begin{eqnarray}&&{\omega }_{j}(t)={z}_{j}(t){\left[{\left(\upsilon \alpha {\mu }^{\alpha }\right)}^{\tfrac{1}{\alpha -1}}\chi (t)+\displaystyle \sum _{i=1}^{{z}_{j}(t)}{\left[{\omega }_{i\to j}(t| {t}_{i})\right]}^{\tfrac{1}{\alpha -1}}\right]}^{\alpha -1}.\end{eqnarray} \tag{ B8 }$$

Notice that equation (B8) cannot be written as equation (B1). Therefore, while the exposures to infectious sources are initially described as isolated events, the agents' memory causes them to become entangled. For instance, for α = 2, equation (B8) can be written as

$$\begin{eqnarray}&&{\omega }_{j}(t)=2\upsilon \mu \chi (t){z}_{j}(t)+{\rm{\Omega }}(t){\left[{z}_{j}(t)\right]}^{2}\end{eqnarray} \tag{ B9 }$$

which has an explicit quadratic dependence on z_j(t).⁴ The simple contagion of the standard SIS model can be recovered only with α = 1, for which equation (B4) equates with equation (B3). In conclusion, the non-Markovianity of the miccSIS model induces an effective social reinforcement/inhibition even though it was not incorporated in the initial description of the model.