A Graph-Based Modelling of Epidemics:
Properties, Simulation, and Continuum Limit

The modelling of infectious disease transmission has a long history in mathematical biology, but in recent years it has gained an increasing influence on the theory and practice of disease management and control. [5, 9, 19, 32, 33]. Most epidemiological models are compartmental models, with the population divided into classes and with assumptions about the rate of transfer from one class to another. In this paper we consider as a basic model the Susceptible-Exposed-Infectious-Removed (SEIR) model describing the disease transmission and the rate of infected individuals. In fact, the SEIR model, and some of its modifications, is widely considered in the literature and in applications as a starting point for describing the spread of various diseases, for example, tuberculosis [46], measles [3], MERS [21], and recently COVID-19 [16]. The basic idea of the SEIR model [32, 33] is to describe the number of infected and recovered individuals based on the number of contacts, probability of disease transmission, incubation period, recovery, and mortality rate. Since the model focuses on a very short time with respect to demographic dynamics, it postulates that births and natural (i.e., not connected with epidemics) deaths balance each other.

While, in some cases, the behaviour of average quantities of a (large) with respect to the time is sufficient to provide useful insight into the spread of the epidemics from the available data, the spatial component of many transmission systems has been recognised to be of pivotal importance. Due to this, spatially heterogeneous features must be included in the model to properly represent the transmission pattern. Then, a reasonable hypothesis about the phenomena may consider that the spatial aspects of transmission heavily influence the aggregation characteristic of the epidemic: we need hence to investigate data by using models that include such spatial connections. For example, the understanding of human mobility and the development of qualitative and quantitative theories is of key importance for the modelling and comprehension of human infectious disease dynamics, on geographical scales of different sizes. Also for the spread of infectious diseases in livestock comprehensive information on livestock movements, cattle movement, and contacts is required to devise appropriate disease control strategies. Understanding contact risk when herds mix extensively, and where different pathogens can be transmitted at different spatial and temporal scales, remains a major challenge [20]. For example, using data related to cattle movements and focusing on the geographical distribution of these movements is possible to improve the analysis of the spread of epizootic diseases [31].

To introduce spatial heterogeneity we consider metapopulation-based models, where the population is partitioned into large, spatially segregated sub-populations. A similar approach could be used in a more general way, irrespective of the biological interpretation: different ages, small interacting communities, and so on [22, 40, 39]. This framework has traditionally provided an attractive approach to incorporating more realistic contact structures into epidemic models since it often preserves analytic tractability but also captures the most salient structural inhomogeneities in contact patterns. Understanding the dynamics of coupled systems on graphs modelling connectivity in real-life systems can be quite challenging. On the other hand, one is often interested in the analysis of average quantities over large networks [36]. For the analysis of the behaviour of the correspondent coupled dynamical systems, we consider the continuum limit where the dynamics on a sequence of large graphs (networks) is replaced by an evolution equation on a continuous spatial domain. Many nontrivial graph sequences have relatively simple limits, described by symmetric measurable functions on a unit square, called graphons. In particular, information on the asymptotic distribution of social contacts can be encoded in the theory of graphons to obtain continuum models that approximate the dynamics over graphs with a large number of nodes.

The paper is organised as follows. In Sect. 2, we briefly previous work on graphons and in Sect. 3, we introduce the SEIR model on a graph and its main properties. Our main goal is to understand the interplay between the transmission of the disease and the distribution of contacts. Our model does not take into account some relevant factors, like the presence of asymptomatic individuals that have a prominent role in disease transmission. We point out that these model improvements could be analysed with methods similar to those introduced in this work. In Sect. 4, we estimate the spectrum of the weighted adjacency matrix involved in the dynamical model and show that the spectrum has a prominent role in the asymptotic behaviour of the epidemic. In Sect. 5, we derive the continuum limit of compartmental epidemiological models on graphs based on a suitable asymptotic procedure, which relies on the notion of graphon. Then, we investigate the relations between the solutions of the graphon-based asymptotics and the SEIR dynamical system on a suitable graph. Finally, we provide some numerical tests. In the numerical experiments and following the detailed analysis in [7], the number of social contacts in the population is described by a Gamma distribution. In Sect. 6, we discuss conclusions and future work.

Graphon theory was introduced and developed in recent years by Lovász, Szegedy, Borgs, Chayes, Sós, and Vesztergombi among others [35, 11, 12, 13, 15]. A graph $G=(V,E)$ is a set of nodes (vertices) $V(G)\,=\,\{1,\,2,\,\ldots,\,n\}$, usually $n=|V(G)|$ where $|\cdot|\leavevmode\nobreak\$ denotes the cardinality of a set, and a set of edges $E(G)\subseteq V\times V$ between the vertices. In what follows the graphs will be simple, without loops or multiple edges, and finite unless otherwise specified. Weights (real numbers) will be given to the edges of a graph to make it an edge-weighted graph. Moreover, we assume that the graph is undirected, i.e., we identify the edges $(i,j)$ and $(j,i)$. Let $\tilde{A}=(\tilde{a}_{ij})$ be the adjacency matrix of the graph:

Let $G$ and $H$ be graphs, a map $\phi$ from $V(H)$ to $V(G)$ is a homomorphism if it preserves edge adjacency, that is if for every edge $(i,j)$ in $E(H)$, $(\phi(i),\,\phi(j))$ is an edge in $E(G)$. Denote by $hom(H,G)$ the number of homomorphisms from $H$ to $G$. For example for $H$ with one node and no edge $hom(H,G)=|V(G)|$, instead when $H$ has two nodes and one edge $hom(H,G)=|E(G)|$, or when $H$ has three vertices and $E(H)=\{(1,2),\,(2,3),\,(3,1)\}$, $hom(H,G)$ is $6$ times the number of triangles in $G$. Normalising by the total number of possible maps, we get the density of homomorphisms from $H$ to $G$, $t(H,G)=\frac{hom(H,G)}{|V(G)|^{|V((H)|}}$. It is defined that a sequence of simple graphs $\{G_{n}\}_{n\in\mathbb{N}}$ is convergent if $G_{n}$ become more and more similar as $n$ goes to infinity.

The breakthrough results of Lovász and Szegedy [35] is that a (left) convergent undirected graph sequence has a limit object, which can be represented as a measurable function.

Let $\cal{W}$ denote the space of all the graphon and $\cal{W}_{0}$ is the set of graphons with values in $[0,1]$. In our models, we will consider, due to the meaning of the connections between different sub-population only graphons in $\cal{W}_{0}$. We give below the representation Theorem of Lovász and Szegedy [34].

Heuristically the intuition behind this definition is that the interval $[0,1]$ represents a sort of continuum of vertices, serving as locations or indices, and $W(x,y)$ denotes the probability of having an edge between vertices $x$ and $y$. It is possible to represent a graph $G$ with a graphon from a suitable family of graphons. If $|V(G)|=n$ and choosing any labelling of the nodes, we divide the interval $[0,1]$ in $n$ sub-intervals $I_{j}^{n}$, $j=1,2,\ldots,n$ as

and define the piecewise constant graphon $W_{G}$ by

Then, $W_{G}$ is a kind of functional representation of the adjacency matrix $A$ of the graph $G$. We point out that different labelling yields different graphon associated with the same graph.

It remains to formalise the convergence of the associated graphons to a graph and understand in what sense the space of graphons completes the space of finite graphs. To do this, we define a metric on the set of labelled graphons. In this work, we will use the set of graphons $\cal{W}_{0}$ but the same definitions can be extended to the general graphon space $\cal{W}$. By analogy with the matrix cut norm, we can define the cut norm of $W\in\cal{W}_{0}$ by setting

where the supremum is taken over all measurable subsets $S$ and $T$. The notion of cut-norm was first introduced by Frieze and Kannan [29], and its fundamental importance for the theory of graphons was recently unveiled in several papers, see in particular [11, 12, 13, 14, 15] and the references therein. In particular, we point out that the notion of cut-norm for the graphons associated with a sequence of graphs is connected with the notion of the (left) convergence for graph sequences, see [14, Theorem 3.8]. Very loosely speaking, a sequence of graphs is left convergent if the local structure of the graphs is asymptotically stable, in the sense that asymptotically the graphs contain the “same number of copies” of any given subgraph. A drawback of the notion of the cut norm is related to the fact that for two graphs $G$, $G^{\prime}$ which are the same but with different node labelling the corresponding adjacency matrices are not the same. Then, we have $\|W_{G}-W_{G^{\prime}}\|_{\Box}\neq 0$. This observation leads to the definition of the cut distance.

Since the cut metric of two different graphons can be zero, strictly speaking, it is not a metric. By identifying graphons $U$ and $W$ for which $\delta_{\Box}(W,U)=0$, we can construct a metric space $(\cal{W}_{0},\delta_{\Box})$ which is compact under this identification [36]. For a sequence of graphons $\{W_{n}\}_{n\in\mathbb{N}}$ in $\cal{W}_{0}$ it is possible to prove [14] that the convergence of $t(H,W_{n})$, as $n\rightarrow\infty$, for any simple graph $H$ is equivalent to requesting that  $\{W_{n}\}_{n\in\mathbb{N}}$ is a Cauchy sequence for $\delta_{\Box}$ or that exists $W\in\cal{W}_{0}$ such that $\delta_{\Box}(W_{n},W)\rightarrow 0$. Using the association between simple graphs and graphons we have the following result [14].

The limit in Theorem 2 is unique up to the identification of graphons with a cut distance equal to zero. Graphons are then the limit objects of converging graph sequences. We will construct and analyse the limit of the dynamical systems on a sequence of simple graphs by using the graphon structure. In particular, we characterise the asymptotic behaviour of the corresponding sequence of SEIR models defined on the simple graphs $G_{n}$. Pioneering works about dynamical systems and graphon were done by G S. Medvedev for the Kurtamoto model [38] and nonlinear heat equation [37]. As compared to [37, 38] we handle the time-dependent case, the biological characteristic of the epidemics and the contact dynamics can change with respect to time. Moreover, we consider not only a scalar dynamical system on the single node but a vectorial one and, from the technical hypothesis, we do not require that $W\in L^{2}([0,1]^{2})$. In recent literature it is possible to find related works using graphons to study epidemics spread, e.g. [49, 50, 51].

When investigating the transmission of infectious diseases, the analysis of the average behaviour of a large population is sufficient to provide useful insight and extract valuable information from the input data. However, the importance of the spatial component of many transmission systems is being increasingly recognised [40]. The main approaches for spatial models concern different scales: an individual-based simulation, a meta-population model, or a network model. Individual-based models explicitly represent every individual and usually assume a variable probability that any infectious host can infect any other susceptible host. Then the model should be able to account for the states of all $N$ individuals in the population in an independent manner, and at the same time, it should allow for arbitrary interactions among them. The analysis of these models is a difficult task the computational cost of numerical simulations is very onerous, and the extraction of the collective behaviours is very complex. Conversely, in the meta-population models the number of individuals at different space locations is in some states. These models often assume that each location is connected to others, with possible variable strengths of connection.

To describe a mathematical model for the spread of infectious disease, one has to make some assumptions about the disease transmission. We consider here, as a basic model, the SEIR compartmental model where individuals are classified into different population groups based on the infection status. The model tracks the number of people in each of the following categories: Susceptibles (individuals that may become infected), Exposed (individuals that have been infected with a pathogen, but due to the pathogen incubation period, are not yet infectious), Infectious (individuals that are infected with a pathogen and may transmit it to others), and Recovered (individual that is either no longer infectious or has been “removed” from the population). Indeed, we consider diseases with a latent phase during which the individual is infected but not yet infectious: a recent example of the application of this type of model is the description of the transmission of the COVID-19 disease [26].

In the scalar case, the total (initial) population, $N$, is categorised into four classes, namely, $s(t)$ (susceptible), $e(t)$ (exposed), $i(t)$ (infected-infectious), and $r(t)$ (recovered), where $t$ is the time variable. The basic assumption for the scalar model is the homogeneously mixing population hypothesis which, roughly speaking, means that a given infectious individual may transmit the disease to any susceptible individual at the same rate. Also, one postulates that all the individuals in the population have the same chances of interacting with each other. People move from S to E based on the number of contacts with I individuals, multiplied by the probability of infection $\beta$, where $\beta i(t)/N$ is the average number of contacts with infection per unit time of one susceptible person. The other processes taking place at time $t$ are: the exposed E become infectious I with a rate $\mu$ and the infectious recover R with a rate $\gamma$. Recovered individuals do not flow back into the S class, as lifelong immunity is postulated. The fractions $1/\mu$ and $1/\gamma$ are the average disease incubation and infectious periods, respectively. We assume that the total population remain constant, i.e.,  $s(t)+e(t)+i(t)+r(t)=N$. Next, we consider the rescaled variables, which for simplicity we do not relabel,  $s(t)/N\rightarrow s(t)$, $e(t)/N\rightarrow e(t)$,  $i(t)/N\rightarrow i(t)$, $r(t)/N\rightarrow r(t)$. The ordinary differential equations (ODEs) governing the SEIR model are then

We consider a meta-population model and a network that is represented by a simple graph $\mathcal{G}=(V,E)$ with $n$ vertices (nodes, regions, patches) and $m$ edges (connections). Each edge is described by a couple of nodes $(u,v)$, $u,v\in V$. We order the nodes and label them with an integer index. We assume that the adjacency matrix of $G$ is an irreducible matrix: there are no isolated and unreachable groups. In node $j$ the corresponding sub–population possesses $N_{j}$ individuals, and $\sum_{j=1}^{n}N_{j}=N$. We assume that individuals can move to a different node, interact with people in that node, and then return to the original one. We denote by $\hat{a}_{jk}$ the total amount of the subpopulation $j$ that “goes” to node $k$ and interacts with the people in that node. We call $\hat{A}$ the matrix with entries $\hat{a}_{jk}$, so that $\sum_{k=1}^{n}\hat{a}_{jk}=N_{j}$, $j=1,\dots,n$. Associated to $\hat{A}$, let $P^{out}$ the probability outgoing matrix with entries $p^{out}_{jk}$, where we denote by $p^{out}_{jk}$ the probability (percentage) that the sub–population $j$ “goes” to node $k$. In addition, we denote by $P^{in}$ the probability incoming matrix with entries $p^{in}_{jk}$, where $p^{in}_{jk}$ is now the probability (percentage) of the sub–population in $k$ that “arrived” from $j$. Finally, let $M_{j}=\sum_{k=1}^{n}\hat{a}_{kj}$ be the total amount of people arrived in node $j=1,\dots,n$, so that $\sum_{j=1}^{n}M_{j}=N$ again. Then, for any $j=1,\dots,n$, $\sum_{k=1}^{n}p^{out}_{jk}=\sum_{k=1}^{n}p^{in}_{jk}=1$. Moreover, we have

where $\operatorname{Diag}(x_{1},\,x_{2},\dots,\,x_{n})$ is the diagonal matrix with the vector $(x_{1},\,x_{2},\dots,\,x_{n})^{T}\in\mathbb{R}^{n}$ on the main diagonal.

Following the (SEIR model) four different discrete classes are considered for the statuses of individuals in each: susceptible, exposed, infectious, and recovered. Let $S_{j}(t),\,E_{j}(t),\,I_{j}(t),\,R_{j}(t)$ the number of individuals in the node $j$ at time $t$, $S_{j}(t)+E_{j}(t)+I_{j}(t)+R_{j}(t)=N_{j}$: we consider a time interval in which we can neglect demographics. Without any interaction with other nodes, within a deterministic approach of the compartmental models, with continuous time $t$, the epidemic dynamics can be described by the system of differential equations in (5):

where the parameter $\lambda$ is the force of infection.

Concerning the behaviour of an epidemic, $\lambda$ is the rate at which susceptible individuals become infected or exposed, and it is a function depending on the number of infectious individuals: it contains information about the interactions between individuals that concur with the infection transmission. If we suppose that the population of $N_{j}$ individuals mix at random, meaning that all pairs of individuals have the same probability of interacting, the force of infection may be computed as:

where $\beta$ is the infectious rate. Then the system state,

We denote by $(s_{j}(t),\,e_{j}(t),\,\imath_{j}(t),\,r_{j}(t))$ the rescaled (percentage) quantities of susceptible, infected, removed at time $t$ at the node $j$ normalised to the number $N_{j}$ of individuals associated with that node. Then, we obtain

where $\dot{\imath}_{j}$ stands for the derivative of the function $\imath_{j}$.

Now, we take a node $j$ that is connected with the other nodes as encoded in matrix $\hat{A}$. Then $S_{j}(t)$ can change due to the contribution of susceptible people that come from $j$ reached an adjacent node $k$ and met infectious people in that node. Then the contribution to $\dot{S}_{j}$ due to the interactions in node $k$ is given by the the $p^{out}_{jk}S_{j}=\hat{a}_{jk}s_{j}$ susceptible people that met a population in node $k$ with a proportion of infectious people given by

Let the vectors ${\bf{X}}(t)=(x_{1}(t),\,x_{2}(t),\,\dots,\,x_{n}(t))^{T}\in\mathbb{R}^{n}$, with $X=S,\,E,\,I,\,R$, the $SEIR$ model on the graph $G$ is the following

where $\hat{B}=P^{out}\operatorname{Diag}(M_{1},\ldots,M_{n})^{-1}{P^{out}}^{T}$, and after rescaling (obtained by a premultiplication with $\operatorname{Diag}(N_{1},N_{2},\ldots,N_{n})^{-1}$),

where $A=P^{out}(P^{(in)})^{T}$.

Owing to the classical Cauchy Lipschitz Picard Lindelöf Theorem, the Cauchy problem obtained by coupling system (9) with the initial data

has a local in-time solution. Since $s_{j},e_{j},i_{j},r_{j}$ represent percentages of individuals, the modelling meaningful range is

for every $t\geq 0$ and every $j=1,\dots,n$. The next result states that the solution of (9) does indeed attain values in the above range.

In the following, we adopt the notation ${\bm{1}}=[1\,1\,\ldots,1]^{\top}$, ${\bm{0}}=[0,\,0,\,\ldots,0]^{\top}$, and, for any vectors $\mathbf{x},\,\mathbf{y}\in\mathbb{R}^{n}$ we write $\mathbf{x}<<\mathbf{y}\,\Leftrightarrow x_{i}<y_{i}\,i=1,2,\ldots,n$, $\mathbf{x}\leq\mathbf{y}\,\Leftrightarrow x_{i}\leq y_{i}\,i=1,2,\ldots,n$ and $\mathbf{x}<\mathbf{y}$ if $\mathbf{x}\leq\mathbf{y}$ but $\mathbf{x}\neq\mathbf{y}$. The asymptotic behaviour of system (9) is described by the following Lemma.

We now focus on the transient behaviour of system (9). If ${\bf s}(t)\geq{\bm{0}}$, then the matrix $B(t)=Diag({\bf s}(t))\,A$ is a non-negative, irreducible matrix. Owing to the Perron-Frobenius Theorem [23], this implies that $B(t)$ has a positive eigenvalue, which we denote by $\lambda_{M}(t)=\rho(B(t))$, equal to the spectral radius. Also, there is an associated positive eigenvector, which we denote by $\mathbf{V}_{M}(t)>>0$. Note that, if ${\bf s}(t)\to 0$, then $\lambda_{M}(t)\to 0$. The next result focuses on the transient behaviour of the solutions of (9).

In this section, we always assume that all the components of ${\bf s}(t)$ are strictly positive. Owing to the proof of Lemma 1, this amounts to assuming that they are strictly positive at $t=0$. The coefficient matrix $B(t):=\operatorname{Diag}({\bf s}(t))A$ of the SEIR model on a general simple graph, can be interpreted as a weighted adjacency matrix, where the weight associated with the edge $(i,j)$ is $b_{ij}(t):=s_{i}(t)a_{ij}$; indeed, the weighted adjacency matrix is not symmetric and the weighted graph is directed. To study the spectrum of $B(t)$ and eventually work with a symmetric weighted adjacency matrix, we introduce its “symmetrisation” (Fig., 1)

where the weight associated with the edge $(i,j)$ is $b^{sym}_{ij}(t):=s_{i}^{1/2}(t)a_{ij}s_{j}^{1/2}(t)$. Furthermore, the matrices $B(t)$ and $B_{sym}(t)$ have the same eigenvalues; in fact,

Let $\mathbf{v}\leavevmode\nobreak\$ be an eigenvector of $B(t)$ with eigenvalue $\nu$. Then,

Indeed, $\mathbf{y}:=\operatorname{Diag}({\bf s}(t))^{-1/2}\mathbf{v}\leavevmode\nobreak\$ is an eigenvector of $B_{sym}(t)$ with eigenvalue $\nu$, i.e., the spectrum of $B(t)$ is real and it uniquely identifies the spectrum of $B_{sym}(t)$, and vice versa. Note that the eigenvectors $\mathbf{y}_{i}:=\operatorname{Diag}({\bf s}(t))^{-1/2}\mathbf{v}_{i}$, $i=1,\ldots,n$ are not orthonormal.

From these properties, we consider the reduced (the vector ${\bf r}$ can be reconstructed from vectors ${\bf s}$, $vli$, ${\bf e}$) SEIR and symmetric SEIR models on a general simple graph (Fig., 2)

Bounds to the spectrum of SEIR on graphs Let $d_{max}$ be the maximum degree of a node in the input graph ${G}$, and let $d_{avg}$ be the average degree of a node in ${G}$. Let $\nu_{1}\geq\nu_{2}\geq\ldots\nu_{n}$ be the eigenvalues of the adjacency matrix $\tilde{A}$ of ${G}$. Then, the bound $d_{avg}\leq\nu_{1}\leq d_{max}$ is applied to estimate the variation of the eigenvalues in linear time. The maximum eigenvalue is computed efficiently through the Arnoldi method [30]. We now derive upper and lower bounds to the eigenvalues of $B(t)$, or equivalently $B_{sym}(t)$ (same eigenvalues):