On bounding exact models of epidemic spread on networks

TL;DR

This paper rigorously proves that the N-Intertwined Mean-Field Approximation (NIMFA) provides an upper bound on the true epidemic spread on networks, using comparison theorems and negative correlation properties.

Contribution

It introduces a rigorous proof that NIMFA overestimates the epidemic process and extends the framework to alternative mean-field closures for weighted and directed networks.

Findings

NIMFA provides an upper estimate of the epidemic process.

The epidemic process exhibits negative correlation.

Results apply to arbitrary weighted and directed networks.

Abstract

In this paper we use comparison theorems from classical ODE theory in order to rigorously show that the N-Intertwined Mean-Field Approximation (NIMFA) model provides an upper estimate on the exact stochastic process. The proof of the results relies on the observation that the epidemic process is negatively correlated (in the sense that the probability of an edge being in the susceptible-infected state is smaller than the product of the probabilities of the nodes being in the susceptible and infected states, respectively), which we also prove rigorously. The results in the paper hold for arbitrary weighted and directed networks. We cast the results in a more general framework where alternative closures, other than that assuming the independence of nodes connected by an edge, are possible and provide a succinct summary of the stability analysis of the resulting more general mean-field models.

Tables

Table 1. Table 1: The relation of the joint and marginal probabilities.

	$⟨I_i⟩$	$⟨S_i⟩$
$⟨I_j⟩$	a	b	q
$⟨S_j⟩$	c	d	1-q
	p	1-p

Equations

$$\dot{⟨I_i⟩}=\tau \sum _{j=1}^N{g}_{ij}⟨S_i{I}_j⟩-\gamma ⟨I_i⟩,$$

$$\dot{⟨S_i{I}_j⟩}$$

$$\dot{⟨I_i{S}_j⟩}$$

$$\dot{⟨I_i{I}_j⟩}$$

$$\dot{⟨S_i{S}_j⟩}$$

$$a+b=q,a+c=p,c+d=1-q,b+d=1-p.$$

$$⟨I_i{I}_j⟩\approx W(⟨I_i⟩,⟨I_j⟩).$$

$$\max (x+y-1,0)\le W(x,y)\le \min (x,y)\text{mbox}forallx,y\in [0,1].$$

$$a=pq,b=(1-p)q,c=p(1-q),d=(1-p)(1-q).$$

$$\dot{X_i}=\tau \sum _{j=1}^N{g}_{ij}(X_j-W(X_i,X_j))-\gamma X_i.$$

$$x_i\le ⟨I_i⟩\le X_i$$

$$⟨S_i{I}_j⟩\le ⟨S_i⟩⟨I_j⟩,⟨I_i{I}_j⟩\ge ⟨I_i⟩⟨I_j⟩,⟨S_i{S}_j⟩\ge ⟨S_i⟩⟨S_j⟩.$$

$$\dot{⟨I_i⟩}\le \tau \sum _{j=1}^N{g}_{ij}⟨S_i⟩⟨I_j⟩-\gamma ⟨I_i⟩.$$

$$\dot{Y_i}=\tau \sum _{j=1}^N{g}_{ij}(1-Y_i)Y_j-\gamma Y_i.$$

$$u\le v,\text{mbox}if{u}_i\le v_i\forall \text{leavevmode}i=1,2,\dots n,u<v,\text{mbox}if{u}_i<v_i\forall \text{leavevmode}i=1,2,\dots n.$$

$$\text{mbox}ifu\le v\text{mbox}and{u}_i=v_i\Rightarrow f_i(u)\le f_i(v)\forall \text{leavevmode}i=1,2,\dots ,n.$$

$${\partial }_j{f}_i\ge 0\forall \text{leavevmode}i,j=1,2,\dots n\text{mbox}andi\ne j.$$

$$\dot{⟨I_i⟩}$$

$$=\tau \sum _{j=1}^N{g}_{ij}(⟨I_j⟩-W(⟨I_i⟩,⟨I_j⟩))-\gamma ⟨I_i⟩+\tau \sum _{N=1}^N{g}_{ij}\left(W(⟨I_i⟩,⟨I_j⟩)-⟨I_i{I}_j⟩\right)$$

$$\le \tau \sum _{j=1}^N{g}_{ij}(⟨I_j⟩-W(⟨I_i⟩,⟨I_j⟩))-\gamma ⟨I_i⟩,$$

$$f_i(x)=\tau \sum _{j=1}^N{g}_{ij}(x_j-W(x_i,x_j))-\gamma x_i.$$

$$\dot{X_i}=\tau \sum _{j=1}^N{g}_{ij}(X_j-\min (X_i,X_j))-\gamma X_i.$$

$$⟨S_i⟩⟨I_j⟩-⟨S_i{I}_j⟩=⟨I_i{I}_j⟩⟨S_i{S}_j⟩-⟨S_i{I}_j⟩⟨I_i{S}_j⟩$$

$$⟨I_i{I}_j⟩+⟨S_i{S}_j⟩+⟨S_i{I}_j⟩+⟨I_i{S}_j⟩=1$$

$$⟨S_i⟩=⟨S_i{I}_j⟩+⟨S_i{S}_j⟩,⟨I_j⟩=⟨S_i{I}_j⟩+⟨I_i{I}_j⟩.$$

$$⟨S_i{I}_j⟩⟨I_i{I}_j⟩+⟨S_i{I}_j⟩⟨S_i{S}_j⟩+⟨S_i{I}_j⟩⟨S_i{I}_j⟩+⟨S_i{I}_j⟩⟨I_i{S}_j⟩=⟨S_i{I}_j⟩$$

$$⟨S_i⟩⟨I_j⟩=⟨S_i{I}_j⟩⟨S_i{I}_j⟩+⟨S_i{I}_j⟩⟨I_i{I}_j⟩+⟨S_i{S}_j⟩⟨S_i{I}_j⟩+⟨S_i{S}_j⟩⟨I_i{I}_j⟩.$$

$$\dot{⟨I_1⟩}$$

$$\dot{⟨I_2⟩}$$

$$\dot{⟨SI⟩}$$

$$\dot{⟨IS⟩}$$

$$\dot{⟨II⟩}$$

$$\dot{⟨SS⟩}$$

$$A(t)=⟨II⟩(t)⟨SS⟩(t)-⟨SI⟩(t)⟨IS⟩(t)$$

$$\dot{A}$$

Full text

1 Institute of Mathematics, Eötvös Loránd University Budapest, Hungary

2 Numerical Analysis and Large Networks Research Group,

Hungarian Academy of Sciences, Hungary

3 School of Mathematical and Physical Sciences, Department of Mathematics, University of Sussex, Falmer, Brighton BN1 9QH, UK

Keywords: mean-field model; exact master equation

AMS classification: 34C12, 34C23, 37G10, 60J28, 90B10, 92D30

$\ast$ corresponding author

email: simonp@cs.elte.hu

1 Introduction

Modelling transmission processes on networks, such as epidemics and rumours, has led to many mathematical challenges [21, 11]. This is mainly due to the high dimensionality of the exact model, which is often a continuous time Markov chain where the size of the state space scales as $m^{N}$, where $m$ is the number of states a node can be in and $N$ is the number of nodes in the network [26, 29]. While theoretically the master equations can be given, their rigorous analysis is out of reach due to the high dimensionality. One approach to deal with this challenge is to consider some ‘clever’ averaging, at node or at population level, and proceed to derive evolution equations for some newly defined average quantities. These however, more often than not, depend on other new average quantities which are of higher order, e.g. the expected number of nodes in a certain state usually depends on the expected number of links/edges in certain states. As a rule of thumb, the dependency between moments is broken by making some closure assumptions where higher order moments are approximated by lower-order ones. This then leads to a low-dimensional system of ordinary differential equations (integro-differential and delay differential equations are also possible) or mean-field model.

The approach above has led to a myriad of mean-field models for $SIS$ (susceptible-infected-susceptible) and $SIR$ (susceptible-infected-recovered) epidemics which are able to capture the average behaviour of epidemics on certain network types (e.g. tree networks and networks build using the configuration model). The most well-known such models are: (a) heterogenous degree-based [22, 6], (b) pairwise [23, 10, 4, 26, 29, 13], (c) effective-degree [16], (d) edge-based compartmental [17, 18], (e) pair-based [24, 25, 12] and (f) N-Intertwined Mean-Field Approximation (NIMFA) [31, 30]. When such models perform well, i.e. output from these agrees well with results from the exact or simulation model, one can proceed to analyse them and to derive analytical results concerning the epidemic threshold and final size, or quasi-equilibrium for $SIS$ epidemics. Such explicit relations between network characteristics and spreading dynamics allow us to better understand how these factors interact and will ultimately lead to more efficient prevention and control measures.

In many cases, mean-field models are validated by simply comparing the results of explicit stochastic network simulation (which stands for the exact model) to output from mean-field models. Such tests are usually performed for a limited number of network types and combination of parameter values. While such heuristic arguments are useful, it is desirable that where possible the difference between the exact and mean-field models is rigorously established using sound mathematical arguments. This endeavour has already led to results proving that under some mild conditions on the degree distribution the edge-based compartmental model is exact for $SIR$ epidemics in the limit of large networks built based on the configuration model [19]. In the case of $SIS$ propagation on a complete graph or on a regular random network the model is a density dependent Markov chain and functional analytic tools can be used to prove that the difference between the output of the mean-field and the exact system scales as $1/N$ for large system size $N$ [5, 14]. Besides estimates on the difference, upper and lower bounds have also been derived for the prevalence obtained from the exact model when $SIS$ or $SIR$ propagation is considered on a complete graph, see [1, 2, 3] These result are valid for graphs with very simple structure, motivating research for finding upper and lower bounds for more complex models.

In this paper, we focus on the NIMFA model for $SIS$ epidemics and we will show that this model provides and upper bound on the exact model on arbitrary weighted and directed networks. In [31] it is claimed that the NIMFA model overestimates the prevalence obtained from the exact system, however, the rigorous proof is not presented there. Here, this is done by using some well-known results from the theory of differential inequalities. The paper is structured as follows. In Section 2 the exact model is formulated and a bottom-up approach is used. This means that the model starts at the level of nodes and focuses on the probability of nodes being either susceptible or infected. Also here, the closure at the level of pairs is generalised beyond simply assuming that the state of nodes at the end of an edge are independent. In Section 3 the main result is presented and it is here where we rigorously prove that NIMFA provides an upper bound on the exact model. In Section 4 we provide a proof, based on dynamical systems arguments, of the fact that epidemics are negatively correlated and point out that this is crucial for the proof of the main result. This is followed by the analysis on the closed model in Section 5. The paper concludes with a short discussion of the main findings and possible extensions.

2 Model formulation

Consider a network with $N$ nodes and assume that no node has an edge to itself, but we allow for node $j$ to have an edge to node $i$ having some weight $g_{ij}$. Typically $g_{ij}=1$ if there is an edge from $j$ to $i$ and [math] otherwise, but the model formulation works for any directed and weighted network and thus we can consider $g_{ij}\in[0,\infty)$ for $i,j=1,2,\dots,N$. We can use the adjacency matrix $G=(g_{ij})_{i,j=1,2,\dots,N}$ to represent the network. We assume that the transmission rate from $j$ to $i$ is $\tau g_{ij}$. The recovery rate at each node is $\gamma$. The probability that node $i$ is infected at time $t$ is denoted by $\langle I_{i}\rangle(t)$. The aim is to derive exact and approximate differential equations for this function.

2.1 Exact model

It can be shown from first principles or from the exact master equations formulated in terms of the probabilities of all $2^{N}$ configurations [24] that $\langle I_{i}\rangle(t)$ satisfies the differential equation

[TABLE]

where $\langle S_{i}I_{j}\rangle(t)$ is the probability that the pair consisting of node $i$ and node $j$ is of type $S-I$ at time $t$. This system is exact but not closed hence further differential equations or a closure is needed to determine the probability $\langle I_{i}\rangle(t)$. The differential equations for the pairs take the following form.

[TABLE]

where $\langle A_{i}B_{j}C_{k}\rangle$ is the probability that the triple consisting of the nodes $i,j$ and $k$ is in the state $A-B-C$, and $(i,j)$ runs over all pairs satisfying $1\leq i<j\leq N$ and in all summations $k$ is different from $i$ and from $j$. Note that for any pair $(i,j)$ we have that the right hand sides of the four differential equations sums to 0, that is the sum $\langle S_{i}I_{j}\rangle+\langle I_{i}S_{j}\rangle+\langle I_{i}I_{j}\rangle+\langle S_{i}S_{j}\rangle$ remains constant in time. If this constant is 1 at the initial time instant, then it remains 1 for all time.

2.2 Closure at the level of pairs

Closure at the level of pairs means that the joint probabilities $\langle S_{i}I_{j}\rangle$, $\langle I_{i}S_{j}\rangle$, $\langle I_{i}I_{j}\rangle$ and $\langle S_{i}S_{j}\rangle$ are expressed in terms of the marginal probabilities $\langle I_{i}\rangle$, $\langle I_{j}\rangle$, $\langle I_{i}\rangle$ and $\langle S_{j}\rangle$. This is always an approximation, as we will show below. However, we wish to rigorously quantify the accuracy of the approximation, or to be able to state whether a model based on closures gives upper or lower bounds for the exact values of the joint probabilities, or indeed if the prevalence from such a closed model is below or above the exact prevalence.

First we rigorously define what a closure relation means. The relation of the joint and marginal probabilities are shown in Table 1, where $a=\langle I_{i}I_{j}\rangle$, $b=\langle S_{i}I_{j}\rangle$, $c=\langle I_{i}S_{j}\rangle$ and $d=\langle S_{i}S_{j}\rangle$. All letters denote probabilities hence $a,b,c,d,p,q\in[0,1]$.

The marginals can be expressed in terms of the joint probabilities as

[TABLE]

It can be immediately seen that $a+b+c+d=1$, hence the four equations are not independent, therefore they cannot be solved for the unknowns $a$, $b$, $c$ and $d$. This shows that the marginals do not determine the joint probabilities. However, once one of them is given then the remaining three are determined by system (6). A closure relation means that one of the joint probabilities is given by a certain algebraic relation involving the marginals, and the others are determined by system (6). Here we define the closure for the $II$ pairs, i.e. $a$ is specified in terms of $p$ and $q$. (One can equivalently express the probabilities of $SI$, $IS$ or $SS$ pairs.) So a closure will be a function $W:[0,1]\times[0,1]\to[0,1]$, for which

[TABLE]

In order to have a solution satisfying $a,b,c,d\in[0,1]$, the function $W$ must satisfy some conditions. On one hand, the first two equations of (6) show that $a\leq p$ and $a\leq q$ must hold, that is $a\leq\min(p,q)$. On the other hand, the third and fourth equations should give a nonnegative value for $d$, hence $1-q-p+a\geq 0$ and $1-p-q+a\geq 0$ must hold, leading to $\max(p+q-1,0)\leq a$. (The $\max$ operation is needed because $p+q-1$ may be negative, when the lower bound for $a$ is zero.) Finally, it is natural to assume that $W$ is a symmetric function, i.e. $W(x,y)=W(y,x)$, leading to the following definition of the closure.

Definition 1

A symmetric function $W:[0,1]\times[0,1]\to[0,1]$ is called a closure relation if

[TABLE]

We note that the most important closure relation is $W(x,y)=xy$ leading to the solution of system (6) in the form

[TABLE]

This corresponds to the case when the states of node $i$ and $j$ are independent, i.e. the joint probability is the product of the marginals. To check that $W(x,y)=xy$ satisfies the inequalities in the definition is an easy exercise and is left to the Reader. One can also immediately see that $W(x,y)=\min(x,y)$ and $W(x,y)=\max(x+y-1,0)$ are proper closure relations.

Once a closure relation is chosen, the closed form of (1) can be given as follows. Since $W(\langle I_{i}\rangle,\langle I_{j}\rangle)$ gives an approximation of $\langle I_{i}I_{j}\rangle$ and $\langle S_{i}I_{j}\rangle+\langle I_{i}I_{j}\rangle=\langle I_{j}\rangle$, one can approximate $\langle S_{i}I_{j}\rangle$ as $\langle I_{j}\rangle-W(\langle I_{i}\rangle,\langle I_{j}\rangle)$. Hence the closed system is

[TABLE]

Solving this system for $X_{i}$ yields an approximation for $\langle I_{i}\rangle$.

The ‘art’ of closing (1) lies in choosing $W$ in a way in which $X_{i}$ is as close as possible to $\langle I_{i}\rangle$. In fact, the only function used in the literature is $W(x,y)=xy$ and the accuracy of the closure has been investigated only numerically. Here our goal is slightly different. We want to find upper and lower bounds to $\langle I_{i}\rangle$, i.e. to introduce closures $w$ and $W$ in such a way that for the corresponding solutions $x_{i}$ and $X_{i}$ of (7) the inequalities

[TABLE]

hold.

We will show that $W(x,y)=xy$ gives an upper bound and $W(x,y)=\min(x,y)$ gives a lower bound. However, improving these bounds by more sophisticated closures remains an open question. First, in the next section, we derive conditions on the closure $W$ ensuring that (7) gives an upper or lower bound.

3 Bounds for the exact system

The derivation of the upper bound is based on the fact that $S-I$ pairs are non-negatively, and $I-I$ and $S-S$ pairs are non-positively correlated [7], that is we have for any $i$ and $j$ that

[TABLE]

We will prove these inequalities in Section 4 and use them to give an upper bound for the solution of system (1). The latter is done below. The first relation leads to the differential inequality

[TABLE]

Based on this inequality let us introduce following system of differential equations, called individual-based model or N-intertwined mean-field approximation (NIMFA) [31, 30].

[TABLE]

Since we have the same right hand side in (9) and (10) with inequality in the first, we expect that the NIMFA approximation yields an upper bound on the exact solution. We will prove this by using the comparison theory of ODEs based on the Kamka-Müller condition which is detailed below.

Consider the ODE $\dot{x}(t)=f(x(t))$ and the differential inequality $\dot{y}(t)\leq f(y(t))$ with a given differentiable function $f:\mathbb{R}^{n}\to\mathbb{R}^{n}$ subject to initial conditions satisfying $y(0)\leq x(0)$. The aim is to find conditions on $f$ ensuring that $y(t)\leq x(t)$ for $t\geq 0$. In what follows, the ordering relation for vectors is used in the following sense:

[TABLE]

A sufficient condition on $f$ for the desired inequality to hold is the Kamke-Müller condition [9, 20], which is equivalent to requiring that

[TABLE]

This condition essentially means that the function in the $i$-th coordinate, $f_{i}$ is increasing in all coordinates $x_{j}$ for $j\neq i$. However, this leads to an alternative sufficient condition, which if satisfied, allows us to call the $f$ function cooperative. More precisely, $f$ is called cooperative if

[TABLE]

It can be shown that if $f$ is cooperative in a convex domain, then it satisfies the Kamke-Müller condition. A detailed and comprehensive study of differential inequalities and comparison theory can be found in the book by Szarski [28]. Cooperative systems generate monotone dynamical systems that are dealt with in the book chapter by Hirsch and Smith [8] and the book by Smith [27]. The main comparison result used here is the following.

Lemma 1

Assume that $f$ satisfies the Kamke-Müller condition. If $\dot{x}(t)=f(x(t))$ and the differential inequality $\dot{y}(t)\leq f(y(t))$ holds and $y(0)\leq x(0)$, then $y(t)\leq x(t)$ for $t\geq 0$.

Using this comparison result we can derive the following theorems about the upper and lower bounds.

Theorem 1

Consider a weighted and directed network $G$, the exact individual-based SIS model given by system (1) and the closed system (7) with a closure $W$ satisfying (besides the conditions in Definition 1)

1. (i)

$y\mapsto y-W(x,y)$* is an increasing function,* 2. (ii)

$\langle I_{i}I_{j}\rangle\geq W(\langle I_{i}\rangle,\langle I_{j}\rangle)$* for all $i,j$.*

Assuming that both models start with identical initial conditions, $\langle I_{i}\rangle(0)=X_{i}(0)$ ($i=1,2,\dots,N$), it follows that $\langle I_{i}\rangle(t)\leq X_{i}(t)$ $\forall\leavevmode\nobreak\ i=1,2,\dots,N$ and $\forall\leavevmode\nobreak\ t\geq 0$.

Proof. We start from the exact system,

[TABLE]

where we have used (ii). The right-hand side of the closed system (7) can be given by the function $f:\mathbb{R}^{N}\to\mathbb{R}^{N}$ with coordinate functions

[TABLE]

It can immediately be seen that the solution of the closed system satisfies the differential equation $\dot{x}(t)=f(x(t))$, and the solution of the exact system satisfies the differential inequality $\dot{y}(t)\leq f(y(t))$, with both systems starting from the same initial condition. According to (i) the coordinate function $f_{i}$ is increasing in the variable $x_{j}$, hence $f$ satisfies the Kamke-Müller condition. Therefore, the general comparison Lemma 1 implies that $\langle I_{i}\rangle(t)\leq X_{i}(t)$ $\forall\leavevmode\nobreak\ i=1,2,\ldots,N$ and $\forall\leavevmode\nobreak\ t\geq 0$.

$\Box$

Now, using the closure $W(x,y)=xy$ we can prove that (10) gives an upper bound. Namely, we have seen that $W(x,y)=xy$ satisfies the conditions in Definition 1. Moreover, $y-xy=y(1-x)$ is increasing in $y$ when $x\in[0,1]$, that is (i) holds. In Section 4 we will prove that (ii) also holds, hence $W(x,y)=xy$ satisfies the conditions of Theorem 1 leading to the following corollary.

Corollary 1

Consider a weighted and directed network $G$, the exact individual-based SIS model given by system (1) and the individual-based closed system (10). Assuming that both models start with identical initial conditions, $\langle I_{i}\rangle(0)=Y_{i}(0)$ ($i=1,2,\dots,N$), it follows that $\langle I_{i}\rangle(t)\leq Y_{i}(t)$ $\forall\leavevmode\nobreak\ i=1,2,\dots,N$ and $\forall\leavevmode\nobreak\ t\geq 0$.

Similarly to Theorem 1, the following result can be proved about the lower bound.

Theorem 2

Consider a weighted and directed network $G$, the exact individual-based SIS model given by system (1) and the closed system (7) with a closure $W$ satisfying (besides the conditions in Definition 1)

1. (i)

$y\mapsto y-W(x,y)$* is an increasing function,* 2. (ii)

$\langle I_{i}I_{j}\rangle\leq W(\langle I_{i}\rangle,\langle I_{j}\rangle)$* for all $i,j$.*

Assuming that both models start with identical initial conditions, $\langle I_{i}\rangle(0)=X_{i}(0)$ ($i=1,2,\dots,N$), it follows that $\langle I_{i}\rangle(t)\geq X_{i}(t)$ $\forall\leavevmode\nobreak\ i=1,2,\dots,N$ and $\forall\leavevmode\nobreak\ t\geq 0$.

It is easy to see that $W(x,y)=\min(x,y)$ satisfies the assumptions of the theorem. Namely, we have seen that $W(x,y)=\min(x,y)$ satisfies the conditions in Definition 1. Moreover, $y-\min(x,y)=\max(y-x,0)=\frac{1}{2}(|y-x|+y-x)$ is increasing in $y$, that is (i) holds. In the previous section we saw that $\langle I_{i}I_{j}\rangle\leq\langle I_{i}\rangle$ and $\langle I_{i}I_{j}\rangle\leq\langle I_{j}\rangle$, hence (ii) also holds. Thus $W(x,y)=\min(x,y)$ satisfies the assumptions of Theorem 2 leading to the following corollary.

Corollary 2

Consider a weighted and directed network $G$, the exact individual-based SIS model given by system (1) and the individual-based closed system

[TABLE]

Assuming that both models start with identical initial conditions, $\langle I_{i}\rangle(0)=X_{i}(0)$ ($i=1,2,\dots,N$), it follows that $\langle I_{i}\rangle(t)\geq X_{i}(t)$ $\forall\leavevmode\nobreak\ i=1,2,\dots,N$ and $\forall\leavevmode\nobreak\ t\geq 0$.

4 $S-I$ pairs are non-negatively correlated

In this section we prove that $S-I$ pairs remain non-negatively correlated in the exact system (1)-(5) if they they are non-negatively correlated initially.

Theorem 3

Let $\langle S_{i}\rangle$, $\langle I_{j}\rangle$ and $\langle S_{i}I_{j}\rangle$ solve system (1)-(5). If $\langle S_{i}\rangle(0)\langle I_{j}\rangle(0)-\langle S_{i}I_{j}\rangle(0)\geq 0$ for all $1\leq i<j\leq N$, then $\langle S_{i}\rangle(t)\langle I_{j}\rangle(t)-\langle S_{i}I_{j}\rangle(t)\geq 0$ holds for all $t>0$ as well.

The proof of the theorem will be divided into several propositions.

Proposition 1

[TABLE]

Proof. We will use the identities

[TABLE]

and

[TABLE]

Using these we obtain

[TABLE]

and

[TABLE]

Taking the difference of the last two equations yields the desired relation.

$\Box$

Theorem 3 will be proved by induction according to the the number of nodes in the network. First, we prove the theorem for a single edge, i.e. for $N=2$, for which system (1)-(5) takes the form

[TABLE]

where we wrote $\langle XY\rangle$ instead of $\langle X_{1}Y_{2}\rangle$ for ease of notation.

Proposition 2

The statement of Theorem 3 holds for a graph with two connected nodes.

Proof. We prove that

[TABLE]

is nonnegative if $A(0)\geq 0$, which proves the statement according to Proposition 1.

Using the differential equations (1)-(5) we can derive a differential equation for the function $A$ as follows.

[TABLE]

where

[TABLE]

Thus $A$ satisfies an inhomogeneous linear differential equation. Multiplying this differential equation with $\exp(-2(\tau+\gamma)t)$ and integrating from [math] to $t$ yields

[TABLE]

The non-negativity of $b(s)$ and $A(0)$ yields that $A(t)\geq 0$ for all nonnegative $t$.

$\Box$

Before proving the theorem in the general case, let us introduce the following notations,

[TABLE]

where we have used Proposition 1, and

[TABLE]

which is the conditional probability that node $i$ is susceptible given that node $k$ is susceptible. The correlation of the conditional probabilities are given as

[TABLE]

where the last term is the conditional probability that node $i$ is susceptible and node $j$ is infected given that node $k$ is susceptible. It is important to note that if $A_{ij}^{k}(t)$ and $A_{ij}(t)$ identical as long as node $k$ is susceptible.

The proof of the theorem is based on the differential equation of $A_{ij}$ that is derived in the proposition below.

Proposition 3

The functions $A_{ij}$ for $1\leq i<j\leq N$ satisfy the following system of inhomogeneous linear differential equations.

[TABLE]

where

[TABLE]

and

[TABLE]

Proof. For the ease of notation, we will write $\langle XY\rangle$ instead of $\langle X_{i}Y_{j}\rangle$ and $\langle I_{k}XY\rangle$ instead of $\langle I_{k}X_{i}Y_{j}\rangle$ and similarly for $\langle X_{i}Y_{j}I_{k}\rangle$, where the indices $i$ and $j$ are fixed throughout the proof.

Differentiating (19) and using the differential equations (1)-(5) yield

[TABLE]

where

[TABLE]

Each term will be simplified separately. Simple algebra leads to

[TABLE]

The expression for $Q_{2}$ can be reduced as follows

[TABLE]

where in the last step the identities $\langle II\rangle+\langle SI\rangle=\langle I_{j}\rangle$ and $\langle II\rangle+\langle IS\rangle=\langle I_{i}\rangle$ were used.

Based on the identities

[TABLE]

the term $Q_{3k}$ can be rewritten as

[TABLE]

The transformations of the term $Q_{3k}$ can be completed by using (20) and (21) as follows.

[TABLE]

where we used Proposition 1 for $A_{ki}$.

Similar transformations lead to

[TABLE]

Substituting the expressions obtained for $Q_{1}$, $Q_{2}$, $Q_{3k}$ and $Q_{4k}$ into (23) we get the desired differential equation for $A_{ij}$.

$\Box$

Now we are ready to prove Theorem 3.

Proof of Theorem 3. We prove that $A_{ij}(t)$ is nonnegative if $A_{ij}(0)\geq 0$, which proves the statement according to Proposition 1. This will be proved by induction according to $N$. The statement for $N=2$ is proved in Proposition 2. Assuming that the statement is true for $N-1$ means that $A_{ij}^{k}\geq 0$, where node $k$ can be regarded as a newly added susceptible node that does not invalidate what we already know. This effectively allows us to extend the model from $N-1$ to $N$ nodes, and thus completing the induction step. Let us now consider system (22) for all $i$ and $j$. This is an inhomogeneous system of linear differential equations in the form

[TABLE]

where $x(t)$ represents the vector containing all $A_{ij}$, $M$ is a cooperative matrix, i.e. its off-diagonal entries are nonnegative, and $g(t)\geq 0$. We prove that $x(0)\geq 0$ (coordinate-wise) implies $x(t)\geq 0$ for all $t\geq 0$. Namely, we have $\dot{x}(t)\geq Mx(t)$, since $g(t)$ is non-negative. Consider the homogeneous system $\dot{y}(t)=My(t)$ with initial condition $y(0)=0$. The solution is obviously $y(t)=0$ for all $t$. On the other hand, this system satisfies the Kamke-Müller condition, hence Lemma 1 yields that $x(t)\geq y(t)$ for all $t\geq 0$, leading to the desired statement.

$\Box$

We note that the statement of the Theorem can be formulated in terms of the $I-I$ pairs. Using that $\langle S_{i}\rangle=1-\langle I_{i}\rangle$ and $\langle S_{i}I_{j}\rangle=\langle I_{j}\rangle-\langle I_{i}I_{j}\rangle$, leads to the following.

Corollary 3

Under the assumptions of Theorem 3 the inequality $\langle I_{i}\rangle(t)\langle I_{j}\rangle(t)\leq\langle I_{i}I_{j}\rangle(t)$ holds for all $t>0$.

5 Analysis of the closed model

The goal in this section is to analyse system (7) from the dynamical system point of view. The closure $W$ satisfies the conditions in Definition 1, moreover, according to Corollary 3 it satisfies $W(x,y)\geq xy$ as well. Thus our aim here is to understand the dynamical behaviour of system (7) when $W$ satisfies

[TABLE]

The lower bound $W(x,y)=xy$ is covered in [15], and this result is presented first.

Theorem 4

Given a directed, weighted, and strongly-connected network, $G$, let $\Lambda_{\text{max}}(G)$ be the largest eigenvalue of the adjacency matrix of $G$. Let the closure in (7) be given as $W(x,y)=xy$. If $\gamma<\tau\Lambda_{\text{max}}(G)$, then a unique endemic (nonzero) steady state exists, and it is stable. Moreover, all of its coordinates are positive. If $\gamma>\tau\Lambda_{\text{max}}(G)$, then there is no endemic steady state and the disease-free steady state is stable.

Concerning the upper bound $W(x,y)=\min(x,y)$, the following result can be easily proved.

Proposition 4

Let $G$ be a directed and weighted network and let the closure in (7) be given as $W(x,y)=\min(x,y)$. Then the only steady state is the disease-free and it is stable.

Thus a transcritical bifurcation occurs when the closure is $W(x,y)=xy$ and there is no bifurcation when $W(x,y)=\min(x,y)$; that is the threshold behaviour disappears when such a crude closure is used. This however raises the question of studying the intermediate regime when $W$ is between the two extremes. Below we give a sufficient condition on closures to ensure that the threshold behaviour is maintained.

We will consider closures where $W$ satisfies

[TABLE]

where $V:[0,1]^{2}\to[0,r]$, with some $r\in(0,1)$, is a continuous function satisfying $V(0,0)=0$ and $xy+V(x,y)\min(x,y)\leq\min(x,y)$. We note that the inequalities in (25) yield a sufficient condition for the existence of the transcritical bifurcation. This means that it may be possible to identify closures that lead to transcritical bifurcation but do not satisfy condition (25).

We note that $W(x,y)=xy$ obviously satisfies this condition with $V(x,y)=0$, and a non-trivial example is $W(x,y)=\sqrt{xy}\min(\sqrt{x},\sqrt{y})$. For the latter, simple calculation shows that there exists a $V(x,y)$ such that this is positive and bounded by a constant $r<1$. Below we will prove that for any choice of $W$ that satisfies condition (25) the same threshold as in Theorem 4 is obtained.

The non-trivial steady state $x\in(0,1]^{N}$ of system (7) is given by

[TABLE]

that will be rearranged using $x_{j}-W(x_{i},x_{j})=x_{j}-x_{i}x_{j}+x_{i}x_{j}-W(x_{i},x_{j})$ as

[TABLE]

where $\alpha=\gamma/\tau$, $(Gx)_{i}$ is the $i$-th coordinate of the vector $Gx$ and

[TABLE]

Expressing $x_{i}$ from (26) we get the fixed point equation $x=T(x)$ for the non-trivial steady state with

[TABLE]

We can immediately see that $T$ maps the unit cube $[0,1]^{N}$ into itself and the origin is its fixed point, representing the disease-free steady state. We will show that in the case $\gamma<\tau\Lambda$, that is $\alpha<\Lambda$, $T$ has a nontrivial fixed point in the interior of the cube, representing an endemic steady state. (Here $\Lambda=\Lambda_{\text{max}}(G)$ is the largest eigenvalue of the adjacency matrix of $G$.) The existence of this fixed point will be verified by using Brouwer’s fixed point theorem on a suitably chosen convex subdomain of the cube not containing the origin. In order to achieve this goal we will need a few auxiliary results.

Proposition 5

For a directed, weighted, and strongly-connected network, given by its adjacency matrix $G$, there exists a positive number $\mu$, for which the following holds. If $(Gx)_{i}<\eta$ and $x_{i}\geq 0$ for all $i=1,2,\ldots N$, then $|x|<\mu\eta$.

Proof. Since the network is strongly connected every column of $G$ contains at least one nonzero entry. Hence $Gx\neq 0$ once $x_{i}\geq 0$ for all $i=1,2,\ldots N$ and $x\neq 0$. Therefore

[TABLE]

If $(Gx)_{i}<\eta$ for all $i=1,2,\ldots N$, then $|Gx|<\eta\sqrt{N}$. On the other hand, $|Gx|=|x|\left|G\frac{x}{|x|}\right|>m|x|$, implying $m|x|<\eta\sqrt{N}$. Hence the statement holds with $\mu=\sqrt{N}/m$.

$\Box$

For the next proposition we introduce a function, $h:[0,+\infty)\to[0,1)$, which is defined by

[TABLE]

Proposition 6

For any $\beta<1/\alpha$ there is a $\omega>0$, such that $h(z)<\omega$ implies $\beta z<h(z)$, when $z\neq 0$.

Proof. One can easily check that $h(0)=0$, $h^{\prime}(0)=1/\alpha$, $h^{\prime}(z)>0$ and $h^{\prime\prime}(z)<0$ for all $z\geq 0$, i.e. $h$ is increasing and concave. Hence a line with slope $\beta<h^{\prime}(0)$ passing through the origin, lies below the graph of $h$ in a sufficiently short interval $(0,z_{0})$. Then the statement holds with $\omega=h(z_{0})$.

$\Box$

Proposition 7

For any $x\in[0,1]^{N}$ we have $T_{i}(x)\geq(1-r)h((Gx)_{i})$.

Proof. According to (25) we have $W(x_{i},x_{j})\leq x_{i}x_{j}+rx_{j}$, hence

[TABLE]

yielding $r(Gx)_{i}\geq F_{i}(x)$, that is $(Gx)_{i}-F_{i}(x)\geq(1-r)(Gx)_{i}$. Therefore (27) leads to

[TABLE]

$\Box$

Proposition 8

Let $\alpha<\Lambda$. Then there is a $\rho>0$, such that $T_{i}(x)<\rho$ for all $i$ implies $(Gx)_{i}<\Lambda T_{i}(x)$, when $(Gx)_{i}\neq 0$.

Proof. Choose a small positive $\varepsilon$, for which $\Lambda(1-\varepsilon)>\alpha$ and introduce $\beta=1/\Lambda(1-\varepsilon)<1/\alpha$. Choose $\omega$ to $\beta$ according to Proposition 6. Choose a positive $\delta$ to $\varepsilon$ according to the continuity of $V$ given in (25), that is $|V(x,y)|<\varepsilon$, when $|x|,|y|<\delta$. Finally, determine $\kappa>0$ from

[TABLE]

where $\mu$ is given in Proposition 5. Now we show that choosing a positive $\rho$ satisfying

[TABLE]

the statement holds.

Using Proposition 7, we get for any $i$ that

[TABLE]

Simple algebra shows that $\kappa>h((Gx)_{i})$ implies $(Gx)_{i}<\alpha\kappa/(1-\kappa)$. Hence, according to Proposition 5, we have $|x|<\mu\alpha\frac{\kappa}{1-\kappa}=\delta$. The choice of $\varepsilon$ ensures that $W(x_{i},x_{j})-x_{i}x_{j}\leq\varepsilon x_{j}$, hence

[TABLE]

Therefore (27) leads to

[TABLE]

Now,

[TABLE]

According to Proposition 6, we get that $h((Gx)_{i})>\beta(Gx)_{i}$, hence

[TABLE]

and this completes the proof.

$\Box$

Now we are ready to prove the existence of the endemic steady state and the presence of a transcritical bifurcation.

Theorem 5

Given a directed, weighted, and strongly-connected network, $G$, let $\Lambda$ be the largest eigenvalue of the adjacency matrix of $G$. Let the closure $W$ in (7) satisfy (25). If $\gamma>\tau\Lambda$ then the origin is the only steady state of the system. In the case $\gamma<\tau\Lambda$, an endemic (nonzero) steady state also exists.

Proof. We first consider the $\gamma>\tau\Lambda$ case and take a steady state $x\in[0,1]^{N}$. According to (26) and using that $F_{i}$ is nonnegative we get

[TABLE]

It is easy to see that for two vectors, $u$ and $v$ with nonnegative coordinates, the inequality $0\leq u_{i}\leq v_{i}$ for all $i$ implies $|u|\leq|v|$. Hence for any nonzero steady state $x\in[0,1]^{N}$ we have

[TABLE]

where we used that $\Lambda$ is the largest eigenvalue of $G$. Hence there is no endemic steady state. We note that this part of the proof only used the fact that $W(x,y)\geq xy$, condition (25) has not been used.

Let us turn to the case $\gamma<\tau\Lambda$. As we mentioned above, we will prove the existence of the endemic steady state by applying Brouwer’s fixed point theorem to the mapping $T$ given in (27) in a suitably chosen domain. The goal is to exclude the origin from this domain, hence we introduce a half-space $S=\{x\in\mathbb{R}^{N}:\ \langle x-a,u\rangle\geq 0\}$ with some appropriately chosen vectors $a,u\in\mathbb{R}^{N}$. Then our domain will be $\Omega=[0,1]^{N}\cap S$. In order to have a nonempty intersection we will choose $a$ from the cube $[0,1]^{N}$ and $u$ will be the unique positive eigenvector of $G$ corresponding to the maximal eigenvalue $\Lambda$, that is $Gu=\lambda u$.

In order to prove the invariance of the domain $\Omega$ it is useful to determine the intersection points of the hyperplane given by $\langle x-a,u\rangle=0$ and the coordinate axes. The intersection point on the $i$-th axis is at $c_{i}=\langle a,u\rangle/u_{i}$. It is easy to see that a point $y\in[0,1]^{N}$ is in $\Omega$ if there is a coordinate $i$, for which $y_{i}\geq c_{i}$. Namely, if $y\in[0,1]^{N}$ is not in $\Omega$, then $\langle y,u\rangle<\langle a,u\rangle=c_{i}u_{i}$ for all $i$, hence $y_{i}u_{i}<c_{i}u_{i}$ leading to $y_{i}<c_{i}$ for all $i$. Now, choose $a\in[0,1]^{N}$ in such a way that for all $i$ we have $c_{i}<\rho$ given in Proposition 8, that is $\langle a,u\rangle/u_{i}<\rho$ for all $i$.

We will prove that $T$ maps $\Omega$ into itself. We know that $T$ maps to $[0,1]^{N}$, hence we only need to prove that $\langle T(x)-a,u\rangle\geq 0$ holds for any $x\in\Omega$. We have previously shown that, if there is an index $i$ for which $T_{i}(x)\geq\rho$, then $T(x)\in\Omega$. Hence we only need to consider the case when $T_{i}(x)<\rho$ for all $i$. In this case we can apply Proposition 8 yielding $T_{i}(x)>(Gx)_{i}/\Lambda$. Hence $T_{i}(x)u_{i}>(Gx)_{i}u_{i}/\Lambda$, leading to

[TABLE]

Thus we proved $\langle T(x)-a,u\rangle\geq 0$, which proves that $T$ maps the convex, compact domain $\Omega$ into itself, hence by Brouwer’s fixed point theorem it has a fixed point in $\Omega$, which is a nontrivial steady state completing the proof of the theorem.

$\Box$

We note that the stability of the steady states was also determined in the special case $W(x,y)=xy$. Further conditions on the smoothness of $W$ would make it possible to generalise the stability result of Theorem 3.8 in [11] for different choices of $W$.

6 Discussion

In Section 4 we proved that the closure $W(x,y)=xy$ satisfies the second assumption of Theorem 1, that is this closure leads to an upper bound of the exact model. Similarly, in Theorem 2 we have shown that $W(x,y)=\min(x,y)$ yields a lower bound of the exact model. However, further work could focus on finding more accurate upper and lower bounds with the possibility of constructing a sequence of closures whose limit is closer in some sense to the exact model.

The analysis of the closed model was presented in Section 5. However, the investigation of the local and global stability for a general closure relation is still an open question. Moreover, we have shown that the qualitative behaviour of the closed system depends on the closure and can be significantly different: it may or may not lead to a transcritical bifurcation. The question of whether a closed system shares the same qualitative features as the exact model is an important one, and ideally the behaviour of the two systems should remain the same. Thus identifying the closure or closures which separate different regimes, those that conserve the qualitative behaviour of the exact system versus those that do not, remains an important direction for further research. One possible step towards achieving this may be to find closures that delimit closed models that have different qualitative behaviour when compared to each other, without considering their relation to the exact model.

Acknowledgements

Péter L. Simon acknowledges support from Hungarian Scientific Research Fund, OTKA, (grant no. 115926).