Disease extinction for susceptible–infected–susceptible models on dynamic graphs and hypergraphs Open Access

We consider a susceptible–infected–susceptible (SIS) disease spread on a dynamically evolving graph, where the dynamic edges represent interactions between individuals in a population over time. We also consider a generalization to the case where a disease spreads on a dynamically evolving hypergraph, where individuals interact in groups of size larger than two. Several previous works have considered SIS spreading on a dynamically evolving graph.^1–12 There, some general conditions were found, which guarantee that the disease will vanish asymptotically. However, in several of those works, these conditions require full exact knowledge of the dynamic graph evaluated at all times. In this work, we motivate and study the case where the dynamic graph is a sequence of i.i.d. random graphs, sampled from a static underlying expected graph, and we show how this assumption allows us to derive simple and practical vanishing conditions. Based on a mean field approximation, we find that, in this setting, the spread of the disease on the underlying expected graph provides an accurate description of the spread of the disease on the dynamically evolving graph. We provide both numerical evidence and theoretical analysis to support this observation. We then extend our investigation to the case of a dynamically evolving hypergraph. There too, we proceed via a mean field approximation, in order to derive a static hypergraph on which the spread of the disease is accurately predicted. We illustrate this approximation through numerical simulations.

The main contributions of this paper can be summarized as follows:

• We compare two different models of an SIS disease spreading on a (fixed or dynamic) graph and identify the one, which is more amenable for extension to the hypergraph case. In particular, we formulate and justify a new model, (18), for the spread of disease over a dynamic hypergraph.

• We propose mean field approximations (14), (17), (20), and (22) to these exact models and illustrate their accuracy via numerical experiments. We also show that these mean field approximations provide spectral conditions, which, if satisfied, guarantee the asymptotic vanishing of the disease in the population—see Theorems V.3 and V.4 and the discussion in Sec. VI.

• We contribute further to the theoretical analysis by estimating the difference between the infection rate of dynamic and static exact and mean field models—see Lemma 5.8.

• In developing the new dynamic hypergraph model, in Sec. III B 1, we set up and analyze a hypergraph analog of the classical Gilbert random graph model.

As mentioned in more detail in Section II B, concepts and results in this area are relevant to many scenarios concerning the spread of information, opinion, or disease; however, to be concrete, we describe the models and analysis in the language of epidemiology.

We first consider a population of $n$ individuals for which we have a nonnegatively weighted, undirected graph represented by a symmetric matrix $W∈Rn×n$ that characterizes pairwise affinity. Here, $0≤wij=wji≤1$ represents the strength of the connection between $i$ and $j$ and, therefore, will be used when we quantify the likelihood of transmission of the infection. In this graph setting, it is natural to have no self-loops; that is, $wii=0$ for $1≤i≤n$⁠, but we note that self-loops will arise when we consider flattening of a hypergraph—see the discussion after Definition 2.3.

The SIS modeling framework places every individual in exactly one of two categories: susceptible or infected. Associated with these categories are a parameter $β>0$⁠, which quantifies the overall virulence of the disease, and a parameter $δ>0$⁠, which quantifies the recovery rate. We will assume throughout that recovery of an infected individual is independent of the state of the system.

In describing various SIS models, to avoid an excess of notation, we will use $pi(t)$ to denote both the exact probability that individual $i$ is infected at time $t$⁠, according to a discrete- or continuous-time process, and an approximation to this quantity arising from, for example, a mean field assumption. By referring to specific models, we hope that the precise meaning of $pi(t)$ is clear at each point in the paper.

In the literature, there are two main approaches to modeling an SIS disease spreading on a graph. One common approach, initiated in Ref. 13, considers the parameters $β$ and $δ$ to be, respectively, the probability that an infected individual infects a neighbor during one unit of time and the probability that an infected individual recovers from the disease during one unit of time. One can then derive a model as follows. Let $h>0$ denote the time between two consecutive steps of the process. The probability $1−pi(t+h)$ that node $i$ is susceptible at time $t+h$ can be expressed as the probability $hδpi(t)$ that it was infected at time $t$ but recovered at time $t+h$⁠, plus the probability $((1−pi(t))∏j=1n(1−hβwij(t)pj(t)))$ that it was susceptible at time $t$ and did not get infected by any neighbor at time $t+h$⁠. This gives the following equation for all $1≤i≤n$⁠:

In Ref. 13, spectral conditions were found, which concern the asymptotic vanishing of the disease. These conditions were derived by linearizing the system, keeping only first order terms, yielding

We could also move to a continuous-time setting by subtracting $pi(t)$⁠, dividing by $h$ and taking the limit $h→0$⁠, yielding

We note that (2) may be viewed as an Euler approximation to (3). Use of (1) or (2) or (3) may be regarded as a mean field approach, where we directly consider the probabilities for the nodes to be infected, rather than keeping track of the individual states of the nodes themselves.

Analyzing the linearized model (2) with $h=1$⁠, it was shown in Ref. 13 that the disease asymptotically vanishes; that is, for all $1≤i≤n$⁠, we have $limt→∞pi(t)=0$ if $ρ(W)<δ/β$⁠, where $ρ(⋅)$ denotes the spectral radius. Hence, $βρ(W)/δ$ may be viewed as a graph-level analog of the classical basic reproduction number.

The authors of Ref. 14 consider an alternative approach based on Markov chains. There, we regard the state of each individual as a stochastic process $(Xi(t))t≥0$ taking values in ${0,1}$⁠, with $0$ for susceptible and $1$ for infected. To this two-state valued process, we associate a rate of the transition matrix, where the transition rates between $0$ and $1$ are given as follows. If $Xi(t)=0$⁠, the rate of infection is given by

Here, in a very similar manner to (2), the overall chance of a new infection is taken to be linearly proportional to the number of currently infected neighbors using the affinity weights. If $Xi(t)=1$⁠, the rate of recovery is given by $δ$⁠. Thus, for every $1≤i≤n$⁠, the stochastic process $(Xi(t))t≥0$ satisfies the Markov property. It is not a Markov process; however, the rates of infection are random variables. In order to use Markov theory to analyze the model further, the authors in Ref. 14 effectuate a mean field approximation by considering instead the expected infection rate $E[λ(Xi(t))]=β∑j=1nwijpj(t)$⁠, where $pi(t):=E[Xi(t)]=P(Xi(t)=1).$ We can now associate with the process $(Xi(t))t≥0$ the transition rate matrix

which makes $(Xi(t))t≥0$ into a continuous-time Markov process. Using Markov theory for continuous-time Markov processes (e.g., Chap. 10 in Ref. 15), we can then deduce that $[pi(t)]$ satisfies the differential equation

Let us note that (5) is very similar to the linearized model in (3), but for the $(1−pi(t))$ factor, which, when multiplied with $∑j=1nwijpj(t)$⁠, induces some second order terms. It is argued in Ref. 14 how the simplification in Ref. 13 of (1) to the linear system (2) does not provide a rigorous mean field approximation of the underlying exact processes ${(Xi(t))t≥0}i=1n$⁠. Nonetheless, we are primarily interested, in this paper, in the derivation of sufficient spectral conditions to guarantee the asymptotic vanishing of the disease. To this end, it is sufficient to consider a linear ODE, which dominates the processes ${(pi(t))t≥0}i=1n$ in order to find sufficient spectral conditions. We note furthermore that the spectral vanishing conditions found in Ref. 13 were identical to those found in Ref. 14.

Another important difference to note between (1) and (5), which is of particular relevance when we wish to capture higher-order interactions between individuals, is that only the Markov chain approach (5) naturally extends to a hypergraph setting. Indeed, the rate of infection $λi(X(t))$ defined in (4) does not restrict us to pairwise interactions. The model in (1) on the other hand, via the product term in the equation, directly uses the assumption that individuals interact pairwise, an assumption that ceases to hold in the hypergraph setting. Thus, previous works, which studied epidemic spreading on hypergraphs (e.g., Refs. 16–20), consider a Markov chain approach, extending model (5). We refer to Refs. 21–23 for further details of recent work on deterministic and stochastic models for dynamics on hypergraphs.

To set up our hypergraph regime, we begin with some definitions.

In particular, if all hyperedges have size two, then the hypergraph $(V,E)$ is a standard graph with nodes $V$ and edge set $E$⁠. From now on, let $n:=|V|$ and $m:=|E|$⁠.

We note that in the graph case, where two nodes can be involved in at most one edge, the off-diagonal elements of the clique expansion matrix $W$ reduce to those of the affinity matrix and the diagonal elements of $W$ contain the node degrees—these diagonal entries may be interpreted as weighted self-loops.

In Ref. 16, the authors introduced a nonlinear function $f$ according to which individuals within a hyperedge may propagate the disease in a nonlinear fashion; here, the likelihood for an individual to become infected is not linearly proportional to the number of infected neighbors. Similar models were further investigated in Refs. 17 and 18. The introduction of the nonlinear function $f$ marks a significant departure from the linearity assumption that is used in typical graph-based models, and it allows us to capture features, such as the majority effect²⁴ in opinion dynamics and the viral load effect²⁵ in disease spreading. We refer to Refs. 16–20 for further justification, noting that there is particular interest in the collective suppression case,^19,20 where the function $f$ is assumed to be concave, and the collective contagion case,^17,18,26 where a virus spreads in a hyperedge only if at least a threshold number of individuals in the hyperedge are infected. Given any such function $f$⁠, we can define the rate of infection of a node $i$ by

In the special case of a graph, the rate in (4) is recovered when $f$ is the identity map. As in the graph setting, we wish to make the rates deterministic in order to apply Markov theory to this process, yielding a mean field approximation model. Two mean field approximations have been proposed and studied in Refs. 17, 19, and 20. In Refs. 19 and 20, spectral conditions were found that guarantee the asymptotic vanishing of an SIS disease spreading on a hypergraph for the different mean field models. One of the insights of Ref. 20 is that the two mean field models, while having slightly different vanishing conditions, provide very similar approximations of the exact model. Hence, we shall only focus on one of the two models, the behavior of the other mean field model being conjectured to provide comparable approximations. We consider the mean field approximation obtained by replacing the infection rate (6) with a deterministic rate involving the processes given in ${(pi(t):=E[Xi(t)]=P(Xi(t)=1))t≥0}i=1n$ to obtain

This ensures that $(Xi(t))t≥0$ is a Markov process, and as in the graph setting, we can apply Markov theory to find an ODE system for ${(pi(t))t≥0}i=1n$⁠. We have

This mean field approximation and the exact model were studied in Theorems 6.5, 8.1, and 9.1 of Ref. 19 where it was found that the disease asymptotically vanishes if

where $cf>0$ is a constant depending on the choice of $f$ [for instance, if $f$ is concave, $cf:=f′(0)$ is a valid choice].

We note that the internal sum in (6) runs from $j=1$ to $n$ and, therefore, includes $j=i$⁠. Removing the $j=i$ term from this sum would not have any effect since the infection rate for node $i$ is only relevant when node $i$ is uninfected; that is, $Xi(t)=0$⁠. However, when we move to mean field models, such as (7), removing the $j=i$ term from the sum would make a difference. In this work, in keeping with previous works, we use the full sum, and we note that computational tests suggest that specifying $j≠i$ produces very similar results.

We are interested in analyzing an SIS disease spread on dynamically evolving graphs and hypergraphs. Several settings for dynamic graphs have been investigated.

In Refs. 1–3, 5–8, 27, and 28, the setting consists of a dynamic graph generated from a finite number of graphs, a typical example being a periodic graph $(W(t))t=0∞$⁠, where there exists $T∈N$ such that for each $k∈N$⁠, $W(k+T)=W(k).$ Then, and in other cases where the dynamic graph is generated by finitely many graphs, a spectral condition for the vanishing of the disease was given as

In Ref. 8, the authors derive spectral conditions for the asymptotic vanishing of an SIS disease spreading on a dynamic graph under various assumptions ensuring that the dynamic graph is generated by finitely many affinity matrices ${A(t)|t∈{1,…,T}}$⁠. They argue that a vanishing condition can be found to be generally related to a weighted average of the spectral radii of the finitely many matrices generating the dynamic graph, the weights being proportional to the frequency of occurrence of each affinity matrix in the dynamic evolution of the graph. For instance, if the graph is periodic, cycling through each affinity matrix one by one following a fixed order, the vanishing condition they derive is given by

We note that the authors in Ref. 8 claim that this sufficient condition can be replaced by the more compact condition

However, since it is not true that the sum of the spectral radii is bounded above by the spectral radius of the sum, but rather that it is bounded below by the spectral radius of the sum when the matrices are Hermitian (the spectral radius then being a norm), which is the case here, we do not see how (11) arises.

From (10), we may derive a simple spectral vanishing condition that was also obtained in Theorem $2$ of Ref. 6,

In general, we note that the vanishing conditions derived in these works require full exact knowledge of the dynamic graph, which is not a realistic assumption in many settings.

Alternatively, one can model a dynamic graph as a family of random graphs indexed by time. This is the approach considered in Refs. 4 and 9–11, and the relationship between discrete- and continuous-time network propagation is discussed in Ref. 29. In Ref. 10, the random edges of the graphs evolve according to a Markov process. There, the authors argue that the spectral vanishing condition provided by the static aggregated graph, for which the affinity matrix can be regarded as the long-time average of the affinity matrices of the random graphs, is less informative about the disease-free state of the process than another condition proved in Theorem $3.10$ of Ref. 10. This may suggest, although this was not tested by the authors of Ref. 10, that the mean field approximation induced by the static aggregated graph does not provide an accurate prediction of the exact dynamic model. This observation is an interesting contrast but is not in contradiction with the results of this paper, where we show that the mean field model induced by the static underlying expected graph (which could be regarded as the analog of the static aggregated graph) provides an accurate approximation of the exact dynamic model in a different dynamic setting (see Secs. III A 2 and V and Fig. 1), and where the associated spectral vanishing condition is sharp.

We define our dynamic graph with a different approach to Refs. 1–3 and 5–8. Similarly to Refs. 4 and 9–11, we suppose that the dynamic graph is generated by random graphs. The main difference in our setting is that we assume independence of the random graphs from one time step to the other. We assume the existence of a static underlying expected graph, from which we sample at each time step, identically and independently, a new random sample matrix to represent the interactions of the individuals of the population. Intuitively, the expected graph contains information about all potential interactions a person may have over that time period (of a time step), each interaction carrying a probability weight proportional to the likelihood of this interaction occurring. To simulate a typical day (or some other relevant time period), we, therefore, draw a random sample using the expected graph, which yields random interactions between individuals based on the expected potential interactions, which could occur. We will refer to this as a snapshot modeling regime for the dynamic interactions.

Let $W¯$ be the affinity matrix of the static underlying expected graph. In practice, in order to build $W¯$ synthetically, we will draw the entries of $W¯$ at random, after which they will be treated as fixed in the model. We assume that $W¯$ is a real symmetric $n×n$ matrix, and for all $1≤i≤n$⁠, $w¯i,j∈[0,1]$⁠. It is reasonable to assume that interactions are sparse in any time period; e.g., we do not expect that a person has an equal probability to visit any place in a city, but rather that they will focus a small number of places, typically in a limited neighborhood. To keep $W¯$ sparse, we use a Gilbert graph model. Recall the Gilbert graph distribution $G(n,p)$⁠, where every (undirected) edge of a random graph on $n$ nodes, has independent probability $p∈[0,1]$ to exist.³⁰ We can view $p$ as a sparsity parameter, tuned as a function of $n$⁠, which allows us to control the expected degree of each node, $(n−1)p$⁠, and to ensure that the generated graph does not have too many connections, which could also make related computations too costly to perform. Let $[xij]i,j=1n∈[0,1]n2$ such that $xij=xji$⁠, and let

The dynamic graph is then generated as an i.i.d. sequence of random graphs, sampled independently and identically from $W¯$ as follows. Letting $T′⊂R+$ be the time domain, for each $t∈T′$⁠, we sample independently $Wij(t):=wij(t)∼Bi(1,wij¯)$ if $i<j$⁠, and we let $wij=wji$ if $i>j$⁠. Informally, we toss a biased coin, with bias determined by $wij¯$⁠, to decide whether the undirected edge is created. We set $wii:=0$ for all $1≤i≤n$⁠. This gives us a dynamic graph $(W(t))t∈T′$⁠. We summarize this construction by writing that each random graph was sampled independently as $W(t)∼Bi(1,W¯)$⁠. In particular, if all the entries of $W¯$ are equal to $1$⁠, then $Bi(1,W¯)=G(n,p)$⁠. To each node $i$⁠, we associate a stochastic process $(Xi(t))t∈T$⁠, where $T′⊂T⊂R+$⁠. In practice, $T′$ is assumed to be discrete, while $T$ may be discrete or continuous but must contain $T′$ as a subset. We can naturally extend the definition of $W$ to the possibly larger time domain $T$⁠, where for all $t∈T$⁠, letting $t0:=max{s∈T′|s≤t}$⁠, we let $W(t):=W(t0).$

Motivated by (4), in this dynamic setting, we consider the infection rate of the process to be

This rate may be approximated as before by

where $pi(t):=E[Xi(t)]=P(Xi(t)=1)$⁠, from which we derive the following ODE system to approximate the expected processes:

In a discrete-time setting, for instance, with constant time increment equal to $1$⁠, we have the corresponding recurrence

This last model, (15), is studied in previous work on dynamic graphs.⁸ In Theorems V.3 and V.4, we show that $E[ρ(W0)]<δ/β$⁠, where $W0∼Bi(1,W¯)$ is a random sample independent from $W(t),$ $t≥0$⁠, is sufficient to ensure that $limt→∞pi(t)=0$⁠, almost surely (a.s.) for all $1≤i≤n$ in the continuous case (14) and the discrete case (15), respectively.

Another approximation of the infection rate (12), more simple than (13), is given by

yielding the following ODE system for the associated mean field approximation in a continuous-time setting:

This approximation has the advantage that it summarizes the spread of the disease using a static graph. Such models are better understood (see Sec. II A), and it is known that if $ρ(W¯)<δ/β$⁠, then the disease vanishes asymptotically. From a computational point of view, running a simulation of the spread of a disease on a dynamically evolving graph can be very costly when we work at the scale of a whole city, with a very high number of nodes and interactions. The mean field approximation (17) suggests that we can obtain an accurate prediction of the disease evolution by considering the spread on a static graph, understood as the underlying expected graph of the model, where weighted edges represent the probabilities for an interaction between two individuals to occur. This mean field model, therefore, provides us with an efficient computational method to simulate the spread of disease when the dynamically changing graph is too large to handle directly.

We compare in Sec. V the spectral vanishing conditions and the infection rates associated with the two mean field models (14) and (17).

Under the assumption that the sparsity parameter satisfies $p=ω(log⁡n/n)$⁠, we show in Lemma 5.8 that for a random sample $W0∼Bi(1,W¯)$⁠, we have $ρ(W0)=ρ(W¯)(1+O(1/pn))$⁠, and we show in Lemma 5.7 that $Sn(i)=E[Sn(i)](1+o(1)),$ where $Sn(i):=β∑j=1wij(t)pj(t)$⁠.

In Sec. IV, we verify numerically that the mean field model (17) induced by $W¯$ provides us indeed with an accurate approximation of the behavior of the exact processes ${(Xi(t))t≥0}i=1n$⁠, thus allowing us to approximate our exact model with a disease spreading on a static graph. We do not show the other mean field approximation (14) in the figures below since the results were visually indistinguishable from those with the more simple mean field model (17).

As discussed in Sec. II, modeling the spread of an SIS disease via a Markov process provides us with a more flexible setting, where we can consider interaction within a hyperedge, i.e., where individuals can interact in groups of size larger than two and contaminate each other in a nonlinear fashion. This more general and realistic setting has already been studied in various works^{16,17,19,20,22,26,31} in the case where the higher-order structure remains static in time.

Here, we investigate instead the setting where the disease spreads on a dynamically evolving hypergraph. This poses an important prior problem: which model should we use to generate a random hypergraph? Indeed, based on the various requirements, one may want to impose on the hypergraph (hyperedges sizes, degree distributions, etc.), various models can be used. Here, we develop what we feel is a natural extension of the Gilbert graph model that we used in Sec. III A 1 to generate the random graphs.

Recalling that $n∈N$ represents the number of nodes (i.e., the number of individuals in the population), we now regard $m∈N$ as the targeted number of hyperedges (i.e., the number of places where individuals can interact). We will assign nodes at random to each available hyperedge.

This definition is similar in construction to the Gilbert graph model. For instance, if $n=m$⁠, the incidence matrix of a random hypergraph in $H(n,m,q)$ corresponds to the affinity matrix of a random directed graph in $G(n,q).$ The motivation to use this model is similar to the intuition behind the model described in Sec. III A at the graph level. We assume that each individual has a list of potential places to visit on any given day, with a certain probability. Drawing a new random incidence matrix independently and identically at each time step simulates a typical day of people visiting places in their neighborhood.

Just as in Sec. III A 1, it is more realistic here to consider that the probabilities for a node to be in a hyperedge (i.e., for an individual to visit a place) vary for every pair $(i,h)∈{1,…,n}×{1,…,m}$ so that visits can focus on a subset of locations. We also wish to keep the connections sparse. We, therefore, let $[xih]i,h=1n,m∈[0,1]nm$⁠, $q$ be a sparsity parameter, and set

We, therefore, consider a weighted $n×m$ incidence matrix $I¯$⁠, which, similarly to $W¯$ being the affinity matrix of the underlying expected graph in Sec. III A, can be thought of as the incidence matrix of the underlying expected hypergraph from which we sample independently and identically an incidence matrix $I(t)$ at every time $t∈T′$⁠.

Our sampling procedure of $(I(t))t∈T′$ from $I¯$ is similar to the sampling procedure of $(W(t))t∈T′$ from $W¯$ in Sec. III A 1. At each time step $t∈T′$ and for each node $i∈{1,…,n}$ and hyperedge $h∈{1,…,m}$⁠, we draw independently $Iih(t)∼Bi(1,I¯ih).$ We can observe that this model does not fix a specific configuration of hyperedges of various sizes, the size of a hyperedge being unknown a priori.

As in the graph case, we then consider the Markov processes ${(Xi(t))t∈T}i=1n$⁠, and we can naturally extend the time domain $T′$ of $I$ to the possibly larger domain $T$ by letting for all $t∈T$⁠, $I(t):=I(t0)$⁠, where $t0:=max{s∈T′|s≤t}.$ Following our discussion of hypergraph models in Sec. II, let $f$ be a nonlinear function representing the manner in which the virus spreads in a hyperedge and define the infection rate of a node $i$ at time $t∈T$ by

Generating the random graphs in Sec. III A 1 as random Gilbert graphs sets all diagonal entries to $0$⁠. On the other hand, note that the sum index $j$ spans through all values in ${1,…,n}$⁠, including $i$ in (18). Assuming that $f=Id$⁠, we recover the infection rate (12) of a dynamic graph, but the induced graph here would be given by $W=IIT$⁠; hence, the diagonal entries would be nonzero, indicating the degree of each node instead. As discussed at the end of Sec. II B, we could, for instance, modify our infection rate (12) by restricting the inner sum to $j∈{1,…,n}∖{i}$⁠, leading to slightly different mean field models. Such a change was found to have negligible impact on the final plots in the numerical experiments in Sec. IV and would not affect the spectral vanishing conditions derived in Sec. V.

Note that $H(n,m,q)$ also induces a probability space on the set of weighted graphs with $n$ nodes. From Definition 2.3, given the incidence matrix of a hypergraph $I$⁠, there corresponds a graph represented by $W:=IIT$⁠, which we could call the weighted clique expansion of the hypergraph. Every entry of this matrix is a random variable defined by

Thus, the entries of the induced affinity matrix $W$ take random values in ${0,1,…,m}$⁠. Little seems to be known about how the choice of the parameters in $H(n,m,q)$ can affect properties of the induced random graph associated with $W$⁠. In practice, it is important to know, for instance, how the sparsity parameter can be tuned to control the expected degree of the nodes and ensure that the graph is connected with high probability. In the case of a random Gilbert graph, it is known that connectivity of the graph holds with high probability if $p=Clog⁡n/n$ for $C>1$⁠.³² One may ask whether knowledge about the well-studied properties of random Gilbert graphs leads to results about the random Gilbert hypergraphs described above and their induced random graph.

Let $W=IIT$⁠, where $I∈H(n,m,q)$⁠, and let $M∈G(n,p)$ be a random Gilbert graph.

By independence, we have

while $E[Mij]=p$⁠. Turning to the variance, we have

while $Var(Mij)=p(1−p).$

Equating the expected values, which is equivalent to equating the expected nodal degree of the two random graphs, we would choose $p:=(q/m)1/2$⁠. Hence, we see that for an equal expected value, the variance of the entries of $W$ is larger than the variance of the entries of $M$ but remains finite and contained in $[0,1]$⁠. We may then argue heuristically that by equating the two expected degrees of the two graph models, the variance does not differ greatly; hence, the connectivity properties of the two graphs should remain similar. We have, therefore, chosen $q:=(Clog⁡n/(mn))1/2$⁠, for some $C>1$⁠, to sparsify the $H(n,m,q)$-induced graph while making sure that it remains connected with high probability. Let us emphasize that this is just a rule of thumb to derive a sensible value for the sparsity parameter $q$ and should not be regarded as a rigorous argument.

The function $f$ in (18) represents the nonlinearity of the infection rate: specifying how the probability for a node to get infected depends on the number of infected group members. A typical choice for $f$ is a concave function; e.g., $x↦log⁡(1+x)$ or $x↦arctan(x)$⁠. Another classic choice is to consider a collective suppression model, where $f$ has the form $x↦c21(x≥c1)$⁠, $c1,c2>0$⁠.

As in Sec. III A, we can consider a mean field approximation with the expected processes ${(pi(t))t∈T}i=1n:=(E[Xi(t)])t∈T$ and the deterministic infection rate at node $i$⁠,

yielding

Following Sec. III A, we propose to further approximate this rate of infection by substituting $Iij(t)$ by its expected value $Iij¯$ for all $t∈T$⁠, thus making the infection rate completely deterministic, yielding at node $i$⁠,

which corresponds to the rate of infection for a disease spreading on a static hypergraph, a setting that was already studied in Refs. 16 and 19, yielding the ODE system

There, it was shown that this model satisfied a similar spectral condition for the vanishing of the disease to the graph-based model. Specifically, if $f′(0)ρ(W¯)<δ/β$⁠, then for all $1≤i≤n$⁠, $limt→∞pi(t)=0,$ where $W¯:=I¯I¯T$⁠. We check the accuracy of this mean field approximation in Sec. IV.

Let us first summarize our approach to simulating the exact individual-level SIS model and its mean field approximation on a dynamic graph. In the simulations, we chose $n=600$⁠. We used i.i.d. weights $wij¯$⁠, where $1≤i<j≤n$⁠, sampled uniformly at random in $[0,1]$⁠, a sparsity parameter $p:=8(log⁡n/n)$⁠, and we let $(aij)i<j$ be i.i.d. random variables where $aij=1$ with probability $p$ and $aij=0$ otherwise. We then built the symmetric matrix $W¯$⁠, where $W¯ij:=wij¯aij$⁠. The sparsity parameter $p$ controls the amount of connectivity between the nodes of the graph. Choosing $p$ too small will make the graph disconnected, while choosing $p$ too large will make the graph behave like a complete graph, increasing the computational cost and failing to capture the levels of interaction that typically arise. We know from the study of random Gilbert graphs³² that the asymptotic connectivity threshold of these random graphs occurs at $log⁡n/n$⁠. Here, we choose $p=8(log⁡n/n)$ so that the connectivity threshold is satisfied even for small $n=600$⁠.

The weights of the matrix $W¯$ indicate the probability for any two individuals $i$ and $j$ in the population to meet at any given time $t$⁠. We let $T′:=N∩[0,T]$⁠, where we choose the final time step to be $T:=200$ so that at every time $t∈T′$⁠, we draw an i.i.d. random sample $W(t)$ from $W¯$⁠, where $W(t)$ is symmetric, and such that for every $i<j$⁠, $Wij(t):=wij(t)=1$ with probability $w¯ij$ and $wij(t)=0$ otherwise.

We then proceed with a standard time discretization of the stochastic processes ${(Xi(t))t∈[0,T]}i=1n$⁠. We fix a small time step $Δt:=0.1$ and advance from time $t$ to time $t+Δt$ at each iteration. As explained in Sec. II C, we extend the definition of $W:=(W(t))t∈T′$ over this finer time partition by setting $W(t):=W(t0)$⁠, where $t0:=max{t′∈T′|t′≤t}.$ We choose a vector $r∈[0,1]n$ of i.i.d. uniformly random values in $[0,1]$⁠, and for every node $i∈{1,…,n}$⁠,

• when $Xi(t)=0$⁠, we set $Xi(t+Δt)=1$ if

  $ri<1−exp⁡(−λi(X(t))Δt)$

  and set $Xi(t+Δt)=0$ otherwise;
• when $Xi(t)=1$⁠, we set $Xi(t+Δt)=0$ if

  $ri<1−exp⁡(−δΔt)$

  and set $Xi(t+Δt)=1$ otherwise.

In Sec. II C, we proposed a mean field model (17) to approximate and simplify the above exact individual-based model. This allows us to consider the spread of a disease on the static underlying expected graph $W¯$ from which we draw the random dynamic graph $W=(W(t))t∈T$⁠, a setting which is well-understood and investigated in previous works (see Sec. II A). We simulated this system of ODE using Euler’s method with time step $Δt=0.1$ and tested numerically how accurately it matches the behavior of the exact process described above.

We simulated the exact individual-level SIS model on a dynamic hypergraph similarly, following Sec. III B. We chose $n=600$ nodes, $m=500$ hyperedges, and following the heuristic in Sec. II B and our choice of $p$ in Sec. IV A, we chose the sparsity parameter as $q:=(8log⁡n/(mn))1/2$⁠. We then created the $n×m$ incidence matrix $I¯$⁠. As in Sec. IV A, we chose i.i.d. weights $x¯ih$ sampled uniformly at random in $[0,1]$⁠. Then, we let $aih:=1$ with probability $q$ and $aih:=0$ otherwise, and we let $Iih¯:=x¯ihaih$ for all $(i,h)∈{1,…,n}×{1,…,m}.$

As in the graph setting, we then discretized the stochastic processes ${(Xi(t))t∈[0,T]}i=1n$ with time step $Δt=0.1.$ The construction of the processes is the same than in Sec. IV A; therefore, we do not repeat it here.

The mean field approximation proposed in Sec. III B follows the ODE system (22) As in Sec. IV A, we simulated this system using Euler’s method with time step $Δt=0.1$⁠.

In these simulations, we fixed the recovery parameter to $δ=1$ and varied the value of the infection parameter $β$⁠, taking all the values in ${(0.02)k|k∈{0,…,15}}$⁠. The plots in Fig. 1 show the proportion of infected individuals after a large time $T=200$ as a function of the infection rate parameter $β$⁠. The blue crosses represent the proportion of infected individuals given by the exact model, averaged over $10$ runs. The red asterisks represent the proportion of infected individuals predicted by our mean field model (17). We observe from Fig. 1 that our mean field model (17) provides a very accurate approximation of the exact model. The vertical line in Fig. 1 intersects the abscissa at $βc:=δ/ρ(W¯)$⁠. In the case where an SIS disease is spreading on a fixed graph with affinity matrix $W¯$⁠, $βc$ represents a threshold value in the sense that if $β<βc$⁠, then the disease vanishes asymptotically from the population (e.g., Refs. 14 and 33). Here, the mean field model is defined on a fixed graph with affinity matrix $W¯$⁠; hence, we can expect that the vertical line in Fig. 1 serves as a threshold value for extinction of the disease, i.e., that the red asterisks will be set to $0$ for $β<βc$ but not necessarily for $β>βc$⁠. On the other hand, there is no such theory available to analyze SIS diseases spreading on dynamic graphs (blue crosses in Fig. 1). We observe nonetheless that due to the close approximation of the exact model by the mean field model (17), $βc$ can also serve as a useful threshold value for extinction for the exact model spreading on the dynamic graph $(W(t))t∈T$⁠.

We also find that the mean field model (22) provides an accurate approximation of the exact individual-based model induced by (18), as shown in Figs. 2–4. Here, as in Sec. IV C 1, we fixed the recovery parameter to $δ=1$ and varied the value of the infection parameter $β$⁠, taking all the values in ${(0.02)k|k∈{0,…,15}}$⁠. The plots show the proportion of infected individuals as a function of $β$ after a large time $T=200$⁠. The blue crosses represent the proportion of infected individuals given by the exact model, averaged over $10$ runs. The red asterisks represent the proportion of infected individuals predicted by our mean field approximation (22). For the nonlinear function $f$⁠, according to which the disease spreads in a hyperedge, we used $f(x):=log⁡(1+x)$ (Fig. 2) and $f(x):=arctan(x)$ (Fig. 3). We also tested the collective contagion model, where we apply $f(x):=x$ for edges, and for each hyperedge of size $k≥3$⁠, we apply $fk(x):=c2(k)1(x≥c1(k))$ (Fig. 4). Here, we chose $c2(k)=c1(k):=(k−1)/2$ for all hyperedges of size $k≥3$⁠. Choosing $c1(k)$ too large, e.g., $c1(k)=k−1$⁠, would prevent the disease from spreading in the hyperedge even for large values of $β$ unless the initial proportion of infected individuals $i0$ is higher than $c1(k)/k$⁠, which would cause the plotted functions in Fig. 4 to remain flat. On the other hand, choosing $c1(k)$ too small would produce a collective contagion model similar to a linear graph-based contagion model. We can observe a discontinuity in the plots shown in Fig. 4, both for the mean field approximation and the exact model, indicating that the predicted long-time proportion of infected individuals in the population is not a continuous function of the infectiousness parameter, $β$⁠. This discontinuity is only observed in the case of collective contagion and contrasts with Figs. 2 and 3, where the nonlinear functions representing the rate of contagiousness in the hyperedges are, respectively, $x↦log⁡(1+x)$ and $x↦arctan⁡(x)$⁠, and where the plots appear to be continuous functions of $β$⁠. In the case of the collective contagion model, the nonlinear functions used to represent the contagiousness of the hyperedges are step functions, with sharply discontinuous changes of value. This implies that a small perturbation in the value of $∑j=1nIjh(t)pj(t)$ or $∑j=1nIjh(t)Xj(t)$ can induce a large change of value in the resulting probability $pj(t)$ for node $j$ to get infected at time $t$⁠. By the construction of the exact model, or by the relation describing the mean field approximation (20), we, therefore, see that a small change in the value of $β$ can recursively induce a large change in the value of the $pj(t)′s$ and, hence, a large change in the final proportion of infected individuals. This effect is the same as the one inducing bistability and hysteresis found in Ref. 17 for the collective contagion model. The effect does not occur in Figs. 2 and 3, however, since the concave functions used in the models are smooth. The spectral vanishing condition in Fig. 4, represented by the green vertical line, is given by $βc:=δ/ρ(W(2)+∑k=3K(c2(k)/c1(k))W(k))$⁠, where $W(k):=I(k)(I(k))T$⁠, with $I(k)$ being the incidence matrix of the sub-hypergraph consisting only of hyperedges of size $k$⁠. We can observe that this condition is less sharp than those in Figs. 2 and 3, which are given by $βc:=δ/f′(0)ρ(W¯)$⁠. A possible explanation for this difference is that $f(x):=log⁡(1+x)$ and $f(x):=arctan⁡(x)$ are concave functions, which can be dominated by a simple linear expression, $f(x)≤f′(0)x$⁠. For the collective contagion analysis, we used the bound $fk(x)≤(c2(k)/c1(k))x$⁠, which replaces discontinuous step functions by smoother concave functions and is likely to introduce greater inaccuracies. We leave for future work the issue of deriving better vanishing conditions for non-concave functions, such as those in the collective contagion model.

In this section, we provide theoretical analysis to support the observations made in our numerical experiments. In previous works on disease spreading on dynamic graphs (Refs. 2, 7, and 8), the analysis focused on the discrete-time mean field approximation (15). We analyze this model and its continuous counterpart (14) with regard to our specific assumptions on the interaction dynamics. In Sec. V A, we derive spectral conditions, which guarantee the asymptotic vanishing of the disease in the population. A corresponding result for the hypergraph case is given in Sec. VI.

In Sec. V B, we provide a comparison between the two mean field approximations (14) and (17), in order to justify why their plots closely match numerically. The latter model has the advantage of considering a static graph $W¯$⁠, instead of a dynamic graph $(W(t))t≥0$⁠, leading to a more straightforward analysis, and as we have seen in Sec. IV, nonetheless provides an accurate prediction of the behavior of the exact model on a dynamic graph.

Here, we derive spectral conditions, which guarantee the asymptotic vanishing of the disease for the mean field approximations (14) and (15). Since we focus on deriving sufficient conditions for the vanishing of the disease, it suffices to find such conditions on a simpler system that dominates the original version. This approach was also used in previous works, as described in Sec. II, where, for instance, the authors of Ref. 13 chose to linearize their ODE system to make it more amenable to analysis. To this end, we appeal to the following result from Ref. 20, which follows as a corollary from Ref. 34.

We also use the following lemma.

Note that there is nothing specific about constructing $W0$ as a matrix of binomial random variables in the above proof; i.e., $W0∼Bi(1,W¯)$⁠. A similar result holds for any other random matrix $W0$ sampled with respect to a fixed expected matrix $W¯$⁠.

A similar argument yields to the same spectral vanishing condition in the discrete-time setting.