Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences
Open AccessResearch articles

Epidemics on hypergraphs: spectral thresholds for extinction

    Abstract

    Epidemic spreading is well understood when a disease propagates around a contact graph. In a stochastic susceptible–infected–susceptible setting, spectral conditions characterize whether the disease vanishes. However, modelling human interactions using a graph is a simplification which only considers pairwise relationships. This does not fully represent the more realistic case where people meet in groups. Hyperedges can be used to record higher order interactions, yielding more faithful and flexible models and allowing for the rate of infection of a node to depend on group size and also to vary as a nonlinear function of the number of infectious neighbours. We discuss different types of contagion models in this hypergraph setting and derive spectral conditions that characterize whether the disease vanishes. We study both the exact individual-level stochastic model and a deterministic mean field ODE approximation. Numerical simulations are provided to illustrate the analysis. We also interpret our results and show how the hypergraph model allows us to distinguish between contributions to infectiousness that (i) are inherent in the nature of the pathogen and (ii) arise from behavioural choices (such as social distancing, increased hygiene and use of masks). This raises the possibility of more accurately quantifying the effect of interventions that are designed to contain the spread of a virus.

    1. Introduction

    Compartmental models for disease propagation have a long and illustrious history [1,2], and they remain a fundamental predictive tool [3,4]. For a stochastic, individual-level model, it has been suggested recently that hyperedge information should be incorporated [57]. Hyperedges allow us to account directly for group interactions of any size, rather than, as in, for example, [810], treating them as a collection of essentially independent pairwise encounters.

    In this work, we contribute to the modelling and analysis of disease spreading on a hypergraph. We include the case where the number of infected nodes in a hyperedge contributes nonlinearly to the overall infection rate; this covers the so-called collective contagion model setting and a new alternative that we call a collective suppression model. The main contributions of our work are

    a mean field approximation ((4.1) and (4.2)) with a spectral condition for local asymptotic stability of the zero-infection state (theorem 6.1) and an extension to global asymptotic stability (theorem 6.5) when the nonlinear infection function is concave,

    for the exact, individual-level model, a spectral condition for exponential decay of the non-extinction probability in the concave case (theorem 8.1) and a spectral bound on the expected time to extinction (corollary 8.2),

    extensions of these results to more general partitioned hypergraph models, where distinct infection rates apply to different categories of hyperedge (theorems 6.2, 6.6 and 8.3 and corollary 8.4),

    results for the non-concave collective contagion model (theorems 9.1 and 9.2),

    a complementary condition that rules out extinction of the disease (theorem 8.6),

    interpretations of these mathematical results: the spectral thresholds for disease extinction naturally distinguish between the inherent biological infectiousness of the disease and behavioural choices of the individuals in the population, allowing us to account for intervention strategies (§10).

    The paper is organized as follows. In §2, we introduce the traditional graph-based susceptible–infected–susceptible (SIS) model and quote a spectral condition that characterizes control of the disease. We then discuss the generalization to hyperedges, and motivate the use of infection rates that do not scale linearly with respect to the number of infected neighbours. Section 3 formalizes the hypergraph model and shows how it may be simulated. In §4, we derive a mean field approximation to characterize the behaviour of the model, and in §5 we give computational results to illustrate its relevance. Section 6 analyses the deterministic mean field setting and gives a spectral condition for long-term decay of the disease. The result is local for general infection rates and global (independent of the initial condition) for the concave case. This spectral condition generalizes a well-known result concerning disease propagation on a graph. Section 7 provides further computational simulations to illustrate the spectral threshold. The full stochastic model is then studied in §8, where we extend the analysis in [8] to our hypergraph setting. Here, we study extinction of the disease in the case where the nonlinearity in the infection rate is concave. We also derive conditions for non-extinction. In §9, we extend our analysis to the so-called collective contagion model proposed in [57]. In §10, we summarize and interpret our results and discuss follow-on work.

    For a review of recent studies of spreading processes on hypergraphs, including the dissemination of rumours, opinions and knowledge, we recommend [11, §7.1.2]. The model that we study fits into the framework of [12]. This work introduced the idea of a nonlinear ‘infection pressure’ from each hyperedge, and derived a mean field approximation that was compared with microscale-level simulation results. In [6], the authors studied this type of model on simplicial complexes of degree up to 2 (a subclass of the more general hypergraph setting) and also studied a mean field approximation. These authors examined the mean field system from a dynamical systems perspective and analysed issues such as bistability, hysteresis and discontinuous transitions. Similarly, in [5,7], a hypergraph version was considered. Our work differs from these studies in (i) focusing on the derivation of spectral thresholds for extinction of a disease in both the exact and mean field settings and (ii) seeking to interpret the results from a mathematical modelling perspective. We mention that it would also be of interest to develop corresponding thresholds for the mean field models in [5,6,12].

    2. Stochastic susceptible–infected–susceptible models

    (a) Stochastic susceptible–infected–susceptible model on a graph

    Classical ODE compartmental models are based on the assumption that any pair of individuals is equally likely to interact—this is the homogeneous mixing case [1]. If, instead, we have knowledge of all possible pairwise interactions between individuals, then this information may be incorporated via a contact graph and used in a stochastic model. Here, each node represents an individual, and an edge between nodes i and j indicates that individuals i and j interact. For a population with n individuals, we may let ARn×n denote the corresponding symmetric adjacency matrix, so nodes i and j interact if and only if Aij=1. In this setting, a stochastic SIS model uses the two-state random variable Xi(t) to represent the status of node i at time t, with Xi(t)=0 for a susceptible node and Xi(t)=1 for an infected node. Each Xi(t) then follows a continuous time Markov process where the infection rate is given by

    βj=1nAijXj(t)2.1
    and the recovery rate is δ. Here, β>0 and δ>0 are parameters governing the strength of the two effects. In this model, we see from (2.1) that the current chance of infection increases linearly in proportion to the current number of infected neighbours.

    This model was studied in [10, theorem 1], where it was argued that the condition

    λ(A)βδ<12.2
    guarantees that the disease will die out. Here, λ(A) denotes the largest eigenvalue of the symmetric matrix A. Further justification for this result may be found, for example, in [8,9]. We note that (2.2) gives an elegant generalization of the homogeneous mixing case (where A corresponds to the complete graph).

    (b) Why use a hypergraph?

    It has been argued [11,1316] that in many network science applications we lose information by recording only pairwise interactions. For example, e-mails can be sent to groups of recipients, scholarly articles may have multiple coauthors and many proteins may interact to form a complex. In such cases, recording the relevant lists of interacting nodes gives a more informative picture than reducing these down to a collection of edges.

    In the setting of an SIS model, we may argue that individuals typically come together in well-defined groups; for example, in a household, a workplace or a social setting. Such groups may be handled by the use of hyperedges, leading to a hypergraph; these concepts are formalized in the next section.

    With a classic graph model, as described in §2a, all types of (pairwise) interactions are treated equally and the rate of infection of a node is linearly proportional to the number of infectious neighbours. With a hypergraph we may consider more intricate contagion mechanisms. For example, using the terminology of [7], the collective contagion model is studied in [57]. Here, infection only starts spreading within a hyperedge after a certain critical number of infectious neighbours has been reached. This type of behaviour is relevant, for example, in a shared office area where infection occurs only when a minimal viral load is exceeded. The critical number may also depend on the environment and group size; in a given office space a small number of workers may be able to socially distance in a way that effectively eliminates the risk of infection. Hence, we may wish to use a different threshold for different sizes and categories of hyperedge.

    We mention here that an alternative type of mechanism may also operate, which we call collective suppression. Imagine that a disease may be contracted through contact with a surface that was previously touched by an infected individual. Now suppose that a group of individuals is likely to use the same physical object, such as a door handle, hand rail, cash machine or water cooler. If an infected individual contaminates the object, then further contamination by other individuals is less relevant. In this case, doubling the number of infected users will increase the risk by a factor less than 2; generally, risk grows sublinearly as a function of the size of the hyperedge.

    These arguments motivate us to allow the rate of infection of a node within a hyperdege to depend on a generic function f of the number of infectious neighbours in that hyperedge; this approach was also taken in [12]. The arguments also motivate us to study the case where the form of the nonlinearity depends on the type of hyperedge. We will be particularly concerned with the case where f is concave, since this is tractable for analysis and allows us to draw conclusions about the collective contagion model. We note that if f is the identity, then we recover linear dependence on the number of infectious neighbours and the hypergraph model is equivalent to a virus spreading on the clique graph of the hypergraph.

    3. Susceptible–infected–susceptible on a hypergraph

    (a) Background

    We continue with some standard definitions [17].

    Definition 3.1.

    A hypergraph is a tuple H:=(V,E) of nodesV and hyperedgesE such that EP(V). Here, P(V) denotes the power set of V.

    We will let n and m denote the number of nodes and hyperedges, respectively; that is, |V|=n and |E|=m. Loosely, a hypergraph generalizes the concept of a graph by allowing an ‘edge’ to be a list of more than two nodes.

    Definition 3.2.

    Consider a hypergraph H:=(V,E). The incidence matrix, I, is the n×m matrix such that Iih=1 if node i belongs to hyperedge h and Iih=0 otherwise.

    It is also useful to introduce W:=IIT. This n×n matrix has the property that Wij records the number of hyperedges containing both nodes i and j. In particular, if H is a graph then W is the affinity matrix of the graph.

    (b) General infection model

    In our context, the nodes represent individuals and a hyperedge records a collection of individuals who are known to interact as a group. As in the graph case introduced in §2a, we use a state vector X(t), which follows a continuous time Markov process, where, for each 1in, Xi(t)=1 if node i is infected at time t and Xi(t)=0 otherwise. We continue to assume that an infectious node becomes susceptible with constant recovery rate δ>0. However, generalizing (2.1), we now assume that a susceptible node i becomes infectious with rate

    βhEIihf(j=1nXj(t)Ijh),3.1
    where β>0 is a constant. In (3.1), f:R+R+ specifies the manner in which the contribution to the overall level of infectiousness from each hyperedge involving node i is assumed to increase in proportion to the number of infected nodes in that hyperedge. Throughout our analysis, we will always assume that f(0)=0 and f is C1 in a neighbourhood of 0.

    If f is the identity, the rate of infection reduces to βj=1nWijXj(t). This gives a weighted version of the infection rate of an SIS model on a graph. As discussed in §2b, it may be appropriate to choose nonlinear f in certain circumstances. We note that in [12] the authors have in mind functions which behave like the identity near the origin and have a horizontal asymptote. Instances of such functions are xarctan(x) and xmin{x,c} for some c>0. Relaxing these conditions, we may ask more generally in such a setting that the function be concave. On the other hand, the authors in [5,6] consider a collective contagion model, where infection spreads within a hyperdge only if a certain threshold of infectious vertices is reached in that hyperedge. A collective contagion model may be represented via the function xc21(xc1) for some c1,c2>0, or xmax{0,xc} for some c>0.

    (c) Partitioned hypergraph model

    Following the discussion in §2b, we also introduce a more general case where we partition the hyperedges into K disjoint categories with each category 1kK having its own distinct rate of infection in response to the number of infected nodes in a hyperedge, represented by a function fk. For example, the categories may correspond to different types of housing, workplaces, hospitality venues or sports facilities, each representing different physical spaces and forms of interaction. We may then represent the infection rate of node i as

    βk=1KhEIih(k)fk(j=1nXj(t)Ijh(k)),3.2
    where we let Iih(k)=1 if i belongs to hyperedge h in the category k and Iih(k)=0 otherwise; so I(k) is the incidence matrix of the subhypergraph consisting of only the hyperedges from category k. We will refer to this as a partitioned hypergraph model.

    In this generalized case, a collective contagion model could be defined by first organizing the hyperedges into categories depending on their size, so that category k is the set of hyperedges of size k+1. A collective contagion model may then be represented, for example, via the functions f1:xx and fk:xc2,k1(xc1,k), k{2,,K}.

    4. Mean field approximation

    A classic approach to studying processes with random infection rates is to develop a mean field approximation for the expected process

    (E[Xi(t)])t0=(P(Xi(t)=1))t0=:(pi(t))t0,
    with deterministic rates. For the model in §3b, Xi(t){0,1} and we have Xi(t)Xi(t)+1 with rate given by (3.1) if and only if Xi(t)=0. Conversely we have Xi(t)Xi(t)1 with rate δ if and only if Xi(t)=1. Making the simplifying assumption that the infection rate for node i is independent of its state, we have
    dpi(t)dt=E[βhEIihf(j=1nXj(t)Ijh)]P(Xi(t)=0)δP(Xi(t)=1).
    Now, in order to obtain a model that just involves pi(t), we introduce a further approximation by interchanging the application of f and the expectation in the first factor on the right-hand side. We arrive at the deterministic mean field ODE
    dP(t)dt=g(P(t)),4.1
    where gi:RnR is defined by
    gi(P(t)):=βhEIihf(j=1npj(t)Ijh)(1pi(t))δpi(t).4.2

    5. Simulations and comparison between exact and mean field models

    Let us emphasize that the approximate infection rates in (4.2) differ in general from the expectation of the random rates in (3.1). When the function f is concave, however, Jensen’s reverse inequality indicates that the rates in (4.2) are greater than the expectation of the rates in (3.1). Hence, in this case the expected quantities pi(t) are overestimated by (4.1) and (4.2). This is fine when we are looking for conditions for the disease to vanish. If f is not concave (e.g. for a collective contagion model) it is less clear a priori how the exact model and mean field ODE compare. In this section, we therefore present results of computational simulations in order to gain insight into the accuracy of our mean field approximation.

    (a) Simulation algorithm

    Before presenting numerical results, we summarize our approach for simulating the individual-level stochastic model, which is based on a standard time discretization; see, for example, [12]. Using a small fixed time step Δt, we advance from time t to t+Δt as follows. First, let r[0,1]n be a random vector of i.i.d. values uniformly sampled from [0,1]. For every node 1in,

    when Xi(t)=0, set Xi(t+Δt)=1 if

    ri<1exp(βhIihf(jXj(t)Ijh)Δt),
    and set Xi(t+Δt)=0 otherwise;

    when Xi(t)=1, set Xi(t+Δt)=0 if

    ri<1exp(δΔt),
    and set Xi(t+Δt)=1 otherwise.

    (b) Computational results

    In the simulations, we chose n=400 nodes with fixed recovery rate δ=1 and a time step of Δt=0.05. We look at results for different choices of (i) the infection strength, β, and (ii) the (independent) initial probability for each node to be infectious, which we denote i0. We simulated the mean field ODE using Euler’s method with time step of 0.05. The largest size of a hyperedge was 5 and we distributed the number of hyperedges for the hypergraph randomly as follows: 300 edges, 200 hyperedges of size 3, 100 hyperedges of size 4 and 50 hyperedges of size 5. To give a feel for the level of fluctuation, the individual-level paths are averaged over 10 runs, each with the same hypergraph connectivity and initial state.

    Figures 13 show results for three concave choices of f; respectively,

    f(x)=min(x,3),

    f(x)=log(1+x),

    f(x)=arctan(x).

    Figure 1.

    Figure 1. Here, f(x)=min(x,3). Red dashed line: mean field approximation from (4.1) to (4.2). Blue solid line: proportion of infected individuals, iXi(t)/n, from the individual-level stochastic model (3.1), averaged over 10 runs.(Online version in colour.)

    Figure 2.

    Figure 2. Here, f(x)=log(1+x). Red dashed line: mean field approximation from (4.1) to (4.2). Blue solid line: proportion of infected individuals, iXi(t)/n, from the individual-level stochastic model (3.1), averaged over 10 runs.(Online version in colour.)

    Figure 3.

    Figure 3. Here, f(x)=arctan(x). Red dashed line: mean field approximation from (4.1) to (4.2). Blue solid line: proportion of infected individuals, iXi(t)/n, from the individual-level stochastic model (3.1), averaged over 10 runs.(Online version in colour.)

    For figure 4, we used a collective contagion model on a partitioned hypergraph. Assigning each hyperedge to a category in {1,2,3,4}, where category k contains the hyperedges of size k+1, we chose the following associated functions to determine the infection rates: f1(x):=x, and for k{2,3,4}, fk(x):=(k1)1(xk1).

    Figure 4.

    Figure 4. Collective contagion model on a partitioned hypergraph. Red dashed line: mean field approximation from (4.1) with (6.2). Blue solid line: proportion of infected individuals, iXi(t)/n, from the individual-level stochastic model (3.1), averaged over 10 runs.(Online version in colour.)

    The four figures show the proportion of infectious individuals as a function of time. In these simulations, and others not reported here, we observe that the initial value i0 does not affect the asymptotic behaviour of the process: the process vanishes or converges to a non-zero equilibrium depending on the value of β but regardless of the value of i0. In figures 13, where f is concave, we know that the mean field model gives an upper bound on the expected proportion of infected individuals in the microscale model. We also see that the mean field model provides a reasonably sharp approximation. Moreover, we see a similar level of sharpness in figure 4 for the collective contagion model, where f is not concave.

    A key advantage of the mean field approximation is that it gives rise to a deterministic autonomous dynamical system for which there exists a rich theory to study the asymptotic stability of equilibrium points. This motivates the analysis in the next section.

    6. Stability analysis

    We provide below spectral conditions which imply that the infection-free solution 0Rn is a locally or globally asymptotically stable equilibrium of (4.1) and (4.2). We will find that local asymptotic stability can be shown with no structural assumptions on f. We will also find that global asymptotic stability follows under the same conditions when f is concave. Our conclusions fit into a framework that generalizes the graph case (2.2): the spectral threshold takes the form

    λ(W)cfβδ<1,
    for some constant cf>0 depending only on the choice of f.

    Throughout this work, to be concrete we let |||| denote the Euclidean norm.

    (a) Local asymptotic stability

    Theorem 6.1.

    If

    λ(W)f(0)βδ<1,6.1
    then0Rnis a locally asymptotically stable equilibrium of (4.1) and (4.2); that is, there exists a positive γ such that||P(0)||<γlimt||P(t)||=0.

    Proof.

    We see that g(0)=0, so 0Rn is an equilibrium for (4.1). It remains to show that this solution is locally asymptotically stable. Appealing to a standard linearization result [18], it suffices to show that every eigenvalue of the Jacobian matrix g(0) has a negative real part. We compute

    gipj0={βhIihIj0hf(jpjIjh)(1pi),j0i,βhIihIj0hf(jpjIjh)(1pi)βhIihf(jpjIjh)δ,j0=i.
    We see that g(0)=βf(0)WδI. This matrix is symmetric and therefore has real eigenvalues. Hence, it suffices that the largest eigenvalue of βf(0)W does not exceed δ, and the result follows. ▪

    Theorem 6.1 extends to the partitioned model in (3.2). In this case, gi(P(t)) in the mean field ODE (4.1) is defined as

    gi(P(t)):=βk=1KhEIih(k)fk(j=1npj(t)Ijh(k))(1pi(t))δpi(t),6.2
    and we let W(k):=I(k)I(k)T.

    Theorem 6.2.

    If

    λ(k=1Kfk(0)W(k))βδ<1,6.3
    then0Rnis a locally asymptotically stable equilibrium of (4.1), withgidefined in (6.2).

    Proof.

    The proof of theorem 6.1 extends straightforwardly. We compute

    gipj0={βkhIih(k)Ij0h(k)fk(jpjIjh(k))(1pi),j0i,βkhIih(k)Ij0h(k)fk(jpjIjh(k))(1pi)βkhIih(k)fk(jpjIjh(k))δ,j0=i,gipj0|P=0={βkWij(k)fk(0),j0i,βkWij(k)fk(0)δ,j0=i,
    and note that
    λ(βkfk(0)W(k)δI)<0λ(kfk(0)W(k))<δβ.
     ▪

    (b) Global asymptotic stability for the concave infection model

    We now show that when f is concave the condition in theorem 6.1 ensures global stability of the zero equilibrium, and hence guarantees that the disease dies out according to the mean field approximation.

    Definition 6.3.

    Given a matrix A, define its symmetric version to be

    A(S):=(A+AT)2.

    Lemma 6.4.

    Suppose that A and B aren×nreal matrices, and suppose that there exists a diagonal matrixΛ such that for all i{1,2,,n},Λii0, and

    A=BΛ.
    Then the largest eigenvalues of A and B satisfyλ(A)λ(B), and the largest eigenvalues ofA(S)andB(S)also satisfyλ(A(S))λ(B(S)).

    Proof.

    Let x be a unit eigenvector associated with λ(A(S)). We have

    2λ(A(S))=xTAx+xTATx=xTBx+xTBTx2xTΛx=i,j=1n(bij+bji)xixj2i=1nΛiixi2xT(B+BT)xmax{yT(B+BT)y|yTy=1}=2λ(B(S)).
    The inequality λ(A)λ(B) may be shown similarly. ▪

    Theorem 6.5.

    Suppose f is concave. If

    λ(W)f(0)βδ<1,
    then0Rnis a globally asymptotically stable equilibrium of (4.1) and (4.2); solimt||P(t)||=0for any valid initial condition(that is, with0p(0)i1).

    Proof.

    From the global asymptotic stability result in [19, lemma 1′], it is sufficient to show that all eigenvalues of the symmetric matrix (g(P))(S) are strictly less than 0, for all P0. We have

    gipj0={βhIihIj0hf(jpjIjh)(1pi),j0i,βhIihIj0hf(jpjIjh)(1pi)βhIihf(jpjIjh)δ,j0=i.

    Letting B denote the n×n matrix given by

    Bij0={βhIihIj0hf(jpjIjh)(1pi),j0i,βhIihIj0hf(jpjIjh)(1pi)δ,j0=i,
    we have g(P)=BΛ, where Λ is the n×n diagonal matrix where, for all i{1,2,}, Λii:=βhIihf(jpjIjh)0. On the one hand, lemma 6.4 now yields λ((g(P))(S))λ(B(S)); on the other hand, note that B+δIg(0)+δI, where we interpret the inequality in a componentwise sense and where we use f(jpjIjh)f(0), since f is concave. Hence B(S)+δIg(0)+δI, and since B(S)+δI has only non-negative entries, appealing to the Perron–Frobenius theorem, we have λ(B(S))λ(g(0)).

    Combining these inequalities and using the spectral condition in the statement of the theorem, we deduce that

    λ((g(P))(S))λ(g(0))=λ(βf(0)WδI)<0,
    as required. ▪

    When f is the identity—which is concave—we are effectively using the standard graph-based model, albeit with weighted edges. The inequality λ(W)β/δ<1 in theorem 6.5 then generalizes the well-known vanishing condition (2.2) found, for instance, in [810]. To our knowledge, previous arguments leading to this vanishing condition for graphs were not completely rigorous. Theorem 6.5 provides a new, rigorous justification for this spectral bound in the case of traditional mean field graph models.

    A straightforward adaptation of the proof of theorem 6.5 yields the following global asymptotic stability result for the more general partitioned model.

    Theorem 6.6.

    Suppose allfkare concave. If

    λ(k=1Kfk(0)W(k))βδ<1,
    then0Rnis a globally asymptotically stable equilibrium of (4.1), withgidefined in (6.2).

    7. Simulations to test the spectral condition

    We now show the results of experiments that test the sharpness of our spectral vanishing condition. Here, we used the concave functions f(x)=2log(1+x) (on figures 5a and 6) and f(x)=arctan(x) (on figures 5a and 7) to construct partitioned models with f1(x)=x and fk(x)=f(x) for all k2. We fixed a hypergraph with n=400 nodes, 400 edges, 200 hyperedges of size 3, 100 hyperedges of size 4 and 50 hyperedges of size 5. At time zero, each node was infected with independent probability i0=0.5 and we used a recovery rate of δ=1. In addition to the mean field ODE, we also simulated the microscale model, averaged over five runs, using the discretization scheme described in §5b, with Δt=0.1.

    Figure 5.

    Figure 5. (a) Infection function based on 2log(1+x). (b) Infection function based on arctan(x). For different choices of infection strength β (horizontal axis), we show the proportion of infected individuals at time t=150 (vertical axis) for the mean field approximation (4.1) with (6.2) in red circles and for the individual-level stochastic model (3.2) in blue crosses. The spectral bound arising from our analysis is shown as a green vertical line.(Online version in colour.)

    Figure 6.

    Figure 6. Results with the logarithmic infection function. (a) Proportion of infected individuals using the mean field model. (b) Proportion of infected individuals using the individual-level model. From bottom to top, the β values used were β{(0.02)k|k{0,,10}}. Cases where β is below the spectral bound are coloured in red.(Online version in colour.)

    Figure 7.

    Figure 7. Results with the arctan infection function. (a) Proportion of infected individuals using the mean field model. (b) Proportion of infected individuals using the individual-level model. From bottom to top, the β values used were β{(0.02)k|k{0,,10}}. Cases where β is below the spectral bound are coloured in red. (Online version in colour.)

    In figure 5, the circles (red) show the corresponding proportion of infected individuals according to the mean field model, ipi(t)/n, at time t=150, for a range of different β between 0 and 0.2: β{(0.01)k|k{0,,20}}. The crosses (blue) show the corresponding proportion of infected individuals from the microscale model, iXi(t)/n, at time t=150. The vertical green line represents the critical value βc:=δ/λ(k=1Kfk(0)W(k)) (we have βc0.0265 on the left and βc0.0490 on the right).

    For the mean field model, we know from theorem 6.6 that β<βc guarantees global stability of the zero-infection state. We see that βc also lies close to the threshold beyond which extinction of the disease is lost in the mean field model. For the individual-level stochastic model, theorem 8.3 below shows that β<βc is also sufficient for eventual extinction of the disease. This is consistent with the results in figure 5.

    In figures 6a and 7a, we show individual trajectories of the proportion of infected individuals, ipi(t)/n, according to the mean field model, for a range of β values. For the same range of β values, the plots on the right of these figures show the corresponding proportion of infected individuals from the microscale model, iXi(t)/n. The curves are coloured in red if the spectral vanishing condition β<βc is satisfied. We see qualitative agreement between the mean field and individual-level models, and extinction for the β values below the spectral threshold.

    Having derived and tested spectral conditions that concern extinction of the disease at the mean field approximation level, in the next section, we study the microscale model directly.

    8. Exact model

    To proceed, we recall our setting where at time zero each node has the same, independent, probability, i0, of being infectious; so P(Xj(0)=1)=i0 for all 1jn. This implies that ni0 is the expected number of infectious individuals at time zero.

    We are interested in the stochastic process iXi(t), which records the number of infected individuals. Our analysis generalizes arguments in [8], which considered a stochastic SIS model on a graph with f as the identity map. We are able to prove results in the case of concave nonlinearites, giving insights into the behaviour of the exact model and also allowing comparison with corresponding results for the mean field approximation. Furthermore, we show in §9 how these results can be extended to the case of collective contagion models.

    (a) Extinction

    Our first result shows that the spectral condition arising from the mean field analysis in theorems 6.1 and 6.5 is also relevant to the probability of extinction in the individual-level model.

    Theorem 8.1.

    Suppose f is concave in the hypergraph infection model (3.1). Then

    P(iXi(t)>0)ni0exp((βf(0)λ(W)δ)t).
    Hence, ifλ(W)f(0)β/δ<1then the disease vanishes at an exponential rate.

    Proof.

    Consider the continuous time Markov process {(Yi(t))t0}i=1n taking values in Nn, with transition of states defined for every 1in and t0 by

    {kk+1,with rateβf(0)jWijYj(t),kk1,with rateδYi(t).
    This new process is introduced here purely for the purpose of analysis. However, it may be interpreted as a disease model where the state of each individual is represented by a non-negative integer that indicates severity of infection. Here, exposure to highly infected individuals raises the chance of an increase in infection severity.

    Suppose also that Xi(0)=Yi(0) for all 1in. Since f is concave,

    βhIihf(jXj(t)Ijh)βhIihf(0)jXj(t)Ijh=βf(0)ijWijXj(t),
    from which we see that Yi stochastically dominates Xi. Hence
    P(iXi(t)>0)P(iYi(t)>0)iqi(t),
    where qi(t):=E[Yi(t)]. In terms of the Chapman–Kolmogorov equation, or chemical master equation [20], we have
    dqi(t)dt=βf(0)jWijqj(t)δqi(t).
    Letting Q(t)=[q1(t),q2(t),,qn(t)]T, this linear ODE system solves to give
    Q(t)=exp(t(βf(0)WδI))Q(0).
    The matrix exp(t(βf(0)WδI)) is symmetric and has spectral radius exp((βf(0)λ(W)δI)t). Hence, in Euclidean norm,
    ||Q(t)||exp((βf(0)λ(W)I)t)||Q(0)||.
    Since ||Q(0)||=ni0 and, by Cauchy–Schwarz,
    iqi(t)n||Q(t)||,
    the proof is complete. ▪

    We deduce, analogously to [8], the following corollary.

    Corollary 8.2.

    Suppose f is concave in the hypergraph infection model (3.1). Let τ denote the time of extinction of the disease and supposeλ(W)f(0)β<δ, then

    E[τ]logn+1δf(0)βλ(W).

    Proof.

    Using theorem 8.1,

    E[τ]=0P(τ>t)dt=0P(iXi(t)>0)dtlognδf(0)βλ(W)+(logn)/(δf(0)βλ(W))nexp((βf(0)λ(W)δ)t)dtlogn+1δf(0)βλ(W).
     ▪

    Likewise, the partitioned case yields the following result.

    Theorem 8.3.

    Suppose everyfkis concave in the partitioned hypergraph model with infection rate (3.2). Then

    P(iXi(t)>0)ni0exp(βλ(k=1Kfk(0)W(k))δ).
    Hence, ifλ(k=1Kfk(0)W(k))β/δ<1then the disease vanishes at an exponential rate.

    We also have the following analogue of corollary 7.2 on the expected time to extinction for the partitioned case.

    Corollary 8.4.

    Suppose everyfkis concave in the partitioned hypergraph model with infection rate (3.2). Let τ denote the time of extinction of the disease and supposeλ(k=1Kfk(0)W(k))β/δ<1, then

    E[τ]logn+1δβλ(k=1Kfk(0)W(k)).

    (b) Conditions that preclude extinction

    So far, we have focused on deriving thresholds that imply extinction. In this subsection, following ideas from [8], we derive a condition under which the disease will persist.

    Note that our analysis does not require the graph associated with W to be connected. The disconnected setting is relevant, for example, when interventions have been imposed in order to limit interactions. We let Δ:=DW denote the Laplacian, and let λc(Δ)>0 denote the smallest non-zero eigenvalue of Δ. We also let emax denote the size of the largest hyperedge.

    Definition 8.5.

    Given a hypergraph H, a function f:R+R+ and a subset of the nodes SV, let

    E(S,f):=iShEIihf(jScIjh),
    where Sc:=VS is the complement of S. Also define for integer 1mn/2
    η(H,m,f):=inf{E(S,f)|S||1|S|m},
    and let η(H,m):=η(H,m,Id).

    Note that when S consists of those nodes for which Xi(t)=0, we can write the infection transition rate of i=1nXi(t) as βE(S,f). More generally, βE(S,f) may be regarded as the rate at which nodes in the set S may be infected by nodes in the remainder of the network. When f=Id and m=n/2, η(H,m) is the Cheeger constant, or isoperimetric number, associated with the weighted graph induced by W=IIT. We may also regard η(H,m,f) as the smallest average infection rate over all subsets consisting of no more than half of the network.

    The next theorem gives a probabilistic lower bound on the time to extinction.

    Theorem 8.6.

    Recall that τ denotes the hitting time of the state 0 for the process(jXj(t))t0in the hypergraph model (3.1). If f is concave andλc(Δ)>(2((emax1)/f(emax1)))(δ/β), then

    P(τ>rm+12m)1re(1+O(rm)),
    wherer:=(emax1)δ/f(emax1)βη(H,m)<1andm:=n/2.

    From standard Cheeger inequalities [21], we know that the Cheeger constant of the graph induced by W satisfies 2η(H,m)λc(Δ); hence, using the assumptions on λc(Δ) in the theorem, we see that r<1 indeed.

    In order to prove this result, we introduce the following lemma.

    Lemma 8.7.

    If f is concave and non-decreasing, then

    f(emax1)emax1η(H,m)η(H,m,f).

    Proof.

    The proof is immediate once we see that, by concavity of f, for all x{0,,emax1}

    f(emax1)emax1xf(x).
     ▪

    Now, to prove theorem 8.6 consider the Markov process (Z(t))t0 valued in {0,,m}, with transition of states given by

    {kk+1with transition ratekβf(emax1)emax1η(H,m),kk1with transition ratekδ.
    This Markov process is stochastically dominated by ((i=1nXi(t))t0), which has the same downward transition rate and an upward transition rate that is at least as large
    kβf(emax1)emax1η(H,m)kβη(H,m,f)βE(S,f),
    where S:={i{1,2,,n}|Xi(t)=0}. Thus, to show theorem 8.6 it suffices to find a suitable lower bound for P(τ^>rm+1/2m), where τ^ is the hitting time of 0 for the process (Z(t))t0. This follows by applying theorem 4.1 of [8] to the process (Z(t))t0.

    9. Collective contagion models

    We now consider the collective contagion models from [57], where infection only starts spreading within a hyperedge once a threshold number of infectious nodes in that hyperedge has been reached. As discussed in §2b, collective contagion models can be represented by nonlinear functions of the form f(x):=max{0,xc} for some c>0, or f(x):=c21(xc1) for some c1,c2>0. In these cases, it is obvious that the zero-infection state for the mean field approximation is locally asymptotically stable (and, indeed, theorem 6.1 applies). However, because the functions are not concave, the theory found in §8 for the exact model does not directly apply. Nonetheless, we can still derive similar spectral conditions for the vanishing of the disease by finding concave functions which serve as upper bounds for f. For instance using c21(xc1)c2c1x1(xc1)+c21(xc1), the bounds in theorem 8.1 and corollary 8.2 lead to the following result.

    Theorem 9.1.

    Suppose thatf(x):=c21(xc1)for somec1,c2>0. Then

    P(iXi(t)>0)ni0exp(βc2c1λ(W)δ).
    In particular, ifλ(W)<(c1/c2)δ/β, then the disease asymptotically vanishes with exponential decay and the extinction time τ satisfies
    E[τ]logn+1δβλ(W)c2/c1.

    Likewise note that max{0,xc}(emax1c)/(emax1)x, where we recall that emax is the largest size of a hyperedge of H. Hence, we deduce the following result.

    Theorem 9.2.

    Suppose thatf(x):=max{0,xc}for somec>0. Then

    P(iXi(t)>0)ni0exp(βemax1cemax1λ(W)δ).
    In particular, if
    λ(W)<(emax1emax1c)δβ,
    then the disease asymptotically vanishes with exponential decay and the extinction time τ satisfies
    E[τ]logn+1δβ((emax1c)/(emax1))λ(W).

    10. Discussion

    In this work, we derived spectral conditions that control the spread of disease in an SIS model on a hypergraph. The conditions have the general form

    βλ(W)cfδ<1,10.1
    where cf>0 is a constant depending on the function f that determines the nonlinear infection rate within a hyperedge.

    We note that in the special case where (i) the hypergraph is an undirected graph and hence W becomes the binary adjacency matrix and (ii) we have linear dependence on the number of infectious neighbours for the infection rate of a node, so f is the identity function, the condition (10.1) reduces to the well-known vanishing spectral condition studied in, for example, [810].

    There are two important points to be made about the general form of (10.1). First, the hypergraph structure appears only via the presence of the symmetric matrix WRn×n. Recall that Wij records the number of times that i and j both appear in the same hyperedge. Such weighted but pairwise information is all that feeds into this spectral threshold. On a positive note, this implies that useful predictions can be made about disease spread on a hypergraph without full knowledge of the types of hyperedge present and the distribution of nodes within them. (For example, when collecting human interaction data it is more reasonable to ask an individual to list each contact and state how many different ways they interact with that contact than to ask an individual to list all hyperedges they take part in.) However, this observation also raises the possibility that more refined analysis might lead to sharper bounds, perhaps at the expense of simplicity and interpretability.

    Our second point is that the new vanishing condition (10.1) neatly separates into three aspects:

    (i)

    The biologically motivated infection parameter, β.

    (ii)

    The interaction structure, captured in λ(W).

    (iii)

    The coefficient cf that arises from modelling the nonlinear infection process. For instance, theorems 6.1 and 6.5 have cf=f(0). In the collective contagion model case f(x)=c21(xc1), theorem 9.1 indicates that we can take cf=c2/c1.

    We may view β as an invariant biological constant that reflects the underlying virulence of the disease and is not affected by human behaviour. The factor λ(W), which arises from the interaction structure, will be determined by regional and cultural issues, including population density, age demographics, typical household sizes and the nature of prevalent commercial and manufacturing activities. Interventions, including full or partial lockdowns, could be modelled through a change in λ(W). The third factor, cf, is strongly dependent upon human behaviour and may be adjusted to reflect individual-based containment strategies such as social distancing, mask wearing or more frequent hand washing.

    This work has focused on modelling, analysis and interpretation at the abstract level, concentrating on the fundamental question of disease extinction. Having developed this theory, it would, of course, now be of great interest to perform practical experiments using realistic interaction and infection data, with the aim of

    calibrating model parameters,

    testing hypotheses about the appropriate functional form of the infection rate,

    testing the predictive power of the modelling framework, especially in comparison with simpler homogeneous mixing and pairwise interaction versions,

    quantifying the effect of different interventions.

    Data accessibility

    Code to run the experiments reported here is available at https://www.maths.ed.ac.uk/dhigham/algfiles.html.

    Authors' contributions

    Both authors contributed to the modelling and analysis, and to the document preparation. H.-L.d.K. implemented the computational tests.

    Competing interests

    The authors have no competing interests.

    Funding

    Both authors were supported by Engineering and Physical Sciences Research Council grant no. EP/P020720/1.

    Footnotes

    Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.

    References