Analytic models for SIR disease spread on random spatial networks

Analytic models for SIR disease spread on random spatial networks Abstract We study the propagation of susceptible-infectious-recovered (SIR) disease on random networks with spatial structure. We consider random spatial networks (RSNs) with heterogeneous degrees, where any two nodes are connected by an edge with a probability that depends on their spatial distance and expected degrees. This is a natural spatial extension of the well-known Erdős–Rényi and Chung–Lu random graphs. In the limit of high node density, we derive partial integro-differential equations governing the spread of SIR disease through an RSN. We demonstrate that these nonlocal evolution equations quantitatively predict stochastic simulations. We analytically show the existence of travelling wave solutions on the RSNs. If the distance kernel governing edge placement in the RSNs decays slower than exponential, the population-scale dynamics are dominated by long-range hops followed by local spread of travelling waves. This provides a theoretical modelling framework for disease propagation on networks with spatial structure, extending analytic SIR modelling capabilities from networks without spatial structure to the important real-life case of networks with spatial structure. Potential applications of this modelling framework include studying the interplay between long-range and short-range transmissions in designing polio vaccination strategies, wavelike spread of the plague in medieval Europe, and recent observations of how epidemics like Ebola evolve in modern connected societies. 1. Introduction The spread of infectious disease through human and animal populations exhibits a range of patterns in space and time. In the pre-vaccination era, Measles in the UK spread out from London in a coherent spatial pattern [1]. Other diseases exhibiting coherent spatial expansion include rabies in raccoons [2, 3], foxes [4, 5] and vampire bats [6] as well as the Black Death of 1347–1351 in Europe, which claimed an estimated 30–50% of the European population [7–9]. These coherent spatial patterns seem the norm in many epizootics as well as historical human epidemics. However, more recent human epidemics exhibit different patterns. The long-range connections in modern societies play increasingly important roles in disease propagation [10]. These lead to the appearance of multiple spreading hotspots of disease. The 2013–2016 epidemic of Ebola in West Africa had significant long-range hops: sequencing of over 1600 Ebola virus genomes reveals a heterogeneous and spatially dissociated collection of transmission clusters of varying size, duration and connectivity [11] (Video 1 of [11] illustrates such a pattern of spread in the Ebola epidemic that we may refer to as ‘hop-and-spread’ dynamics). Similarly, SARS and pandemic 2009 H1N1 influenza showed significant local spread, but the global dynamics were dominated by sporadic long-range transmission events. A further example is the influence of spatial structure on polio eradication. Oral polio vaccine (OPV) is a live, weakly transmitted virus which—such as the wild virus—can cause paralysis and can evolve to be more transmissible [12]. Wild type-2 polio is extinct, so ideally we would like to remove OPV-2 vaccination, but cessation comes with the risk that a transmission chain of the vaccine strain may persist until enough unvaccinated children are born to allow it to cause long-scale epidemics [13]. Routine OPV-2 vaccination has stopped, and the vaccine is currently being used in response to identified circulation of vaccine-derived strains. This may eliminate the circulating strain, but with a risk of seeding new transmission chains. Understanding the interplay between long-range and short-range transmissions is important in designing strategies for localized OPV-2 use, and in informing decisions about stockpile maintenance. Network models are a valuable tool for investigating how a population’s contact structure governs the spread of infectious disease [14–16]. They also appear in many other contexts, such as understanding activation patterns in neurons. In many networks of interest, there is an underlying geometry to the location of nodes, and nearby nodes are more likely to be connected than nodes that are separated by a large distance [17–21]. Outside of first-passage percolation [22] where the networks are generally lattices, there is relatively little theoretical study of how a network’s spatial structure affects spreading processes. Commonly used spatial network models have a number of weaknesses. In particular, they are not amenable to analytic investigation, and many of them have subtle preferences for specific directions. By way of contrast, non-spatial network models such as Chung and Lu [23] and Configuration model networks [24] provide significant insight. This is largely because their structure has allowed for analytic exploration. They are ‘locally tree-like’: for fixed cycle length $$D$$, the probability that a random node is in a cycle of length $$D$$ or less goes to zero as the number of nodes grows. The dynamics of many spreading processes, for example, complex contagions [25], the generalized epidemic process [26] and the Watts threshold model [27] can be studied exactly in locally tree-like networks [28], and even when the tree-like hypothesis breaks down, spreading processes on networks are often accurately predicted by the unclustered assumptions [29]. To study these in networks with spatial structure, we would like a spatial network model that contains the locally tree-like behaviour as a limiting case. The susceptible-infectious-recovered (SIR) model is typically formulated as a system of ordinary differential equations (ODEs) for well-mixed populations [30, 31]. To understand how network effects such as social heterogeneity and nonzero partnership duration affect spreading dynamics, the SIR model is often simulated stochastically [32–34]. Recently, [35] introduced the edge-based compartmental modelling (EBCM) technique that leads to low-dimensional ODE systems that match the stochastic dynamics in the limit of large network size for random configuration model networks [24] (also known as Molloy–Reed graphs [36]). It has been shown mathematically in a rigorous way that the ODE approximation of [35] is indeed exact in the limit [37, 38]. The ODE models allow analytic investigation of the SIR disease processes on networks. The networks of [35], however, do not have spatial structure, and there are many important aspects of SIR disease propagation related to the geographical spread of infection that depend on the spatial aspects of population structure [39–41], in addition to the network structure. However, to date, studies of the spread of disease in spatial networks have relied on simulation because the network models used do not lead to a simple analytic model. In this article, we fill this gap, and extend the analytic edge-based compartmental modelling capabilities of [28, 35] for non-spatial networks to spreading processes on networks with spatial structure. We use a spatial network model that forms a natural extension of the well-known Erdős–Rényi and Chung–Lu random graphs, adding the crucial element of spatial structure while keeping the model amenable to mathematical analysis. This allows us to analytically study how spatial structure affects spreading dynamics. In particular, it enables studying the interplay between long-range and short-range transmissions of SIR disease on networks with spatial structure, and how this affects travelling waves and hop-and-spread dynamics, with applications to understanding propagation of diseases like the plague in medieval Europe and Ebola in modern connected societies, and designing polio vaccination strategies. The remainder of this article is structured as follows. In Section 2, we investigate analytic models for SIR disease propagation on random spatial networks (RSNs). We discuss the RSN class of random networks with spatial structure that we consider in this article, and derive an analytic partial integro-differential equation (PIDE) model to quantitatively describe SIR dynamics on RSNs in the limit of large node density. We show analytically that SIR disease propagation on RSNs in this limit allows for exact travelling wave solutions (i.e., nonlinear waves that retain their shape while propagating) characterised by the nonlocal evolution equations of PIDE type. We study the analytic properties of these travelling wave solutions and conditions under which they exist. In Section 3, we apply stochastic SIR simulation to examples of RSNs, and demonstrate that the nonlocal evolution equations quantitatively predict stochastic simulations. We identify travelling waves in the simulations and verify the analytic properties derived from the PIDE. The interplay between short-range and long-range connections reveals dynamics in which travelling waves emanate from focal points that result from long-range hops. Section 4 discusses the application of these new models to problems in disease propagation and concludes. 2. SIR disease propagation on RSNs: analytic models and travelling waves We first introduce the class of RSNs that we consider in this article. We then derive differential equations that govern the spatial dynamics of SIR disease on RSNs in the limit of large node density and demonstrate the existence of nonlinear travelling wave solutions that retain their shapes. We derive analytic properties of the resulting differential equations, including properties of the travelling waves. These models are used in Section 3 to study the dynamics of SIR disease on RSNs. 2.1 The RSN class To study disease propagation on spatial networks, we consider a class of RSNs modelled after Chung–Lu networks [23] (formally a special case of the inhomogeneous random graph class [42]). We place nodes uniformly at random into some region $$V$$ with density $$\rho$$. For our purposes $$\rho = N/|V|$$, where $$N$$ is the total number of nodes and $$|V|$$ is the (assumed finite) volume of the space $$V$$. Each node $$u$$ is independently assigned an expected degree $$\kappa_u$$ from some distribution $$P(\kappa)$$. The average of $$\kappa$$ is denoted $${\left\langle {\kappa} \right \rangle}$$. A $$u$$–$$v$$ edge is created with probability   \begin{equation} p_{uv} = \min\left(\kappa_u\kappa_v \frac{f(d_{uv})}{\rho {\left\langle {\kappa} \right \rangle}}, 1\right). \label{eq:p} \end{equation} (1) We assume that the distance kernel$$f$$ depends only on distance, is non-negative, and integrates to $$1$$: $$\int_V f(|\mathbf{x}-\mathbf{x}_0|) \, \mathrm{d}\mathbf{x} = 1$$ for all $$\mathbf{x}_0$$. Typically, we also assume that $$f$$ decreases monotonically with distance. An alternate function that would generally yield similar behaviour in the large $$\rho$$ limit would be $$p_{uv} = 1- e^{-\kappa_u \kappa_v f(d_{uv})/\rho {\left\langle {\kappa} \right \rangle}}$$, following [43]. If $$p_{uv}<1$$ always then the expected degree of a node $$u$$ is   \[ \int_V \int_0^\infty \rho f(|\mathbf{x}_v-\mathbf{x}_u|) \kappa_u \kappa_v \frac{P(\kappa_v)}{\rho {\left\langle {\kappa} \right \rangle}} \mathrm{d}\kappa_v \, \mathrm{d}\mathbf{x}_v = \kappa_u \, . \] In the large network limit, the degree of $$u$$ is Poisson distributed with mean $$\kappa_u$$. This model offers considerable flexibility: The distance kernel $$f$$ can be tuned to generate a wide range of spatial structure. The expected degree distribution allows us to tune the distribution of realized degrees. Many further generalizations (presented in the discussion in Section 4) emerge naturally. Figure 1 shows part of an RSN with $$10^4$$ nodes in $$[0,1]\times[0,1]$$ with an imposed bimodal degree distribution and a Gaussian distance kernel. Fig. 1. View largeDownload slide An example RSN and its properties. The distance kernel is a Gaussian, $$f(d) = \exp(-d^2/2\sigma^2)/2\pi\sigma^2$$ with $$\sigma = 0.03$$. The imposed distribution of expected degrees is $$P(2)=P(15)=0.5$$. The density is $$\rho=10000$$. One node and its neighbours are highlighted. A random network without spatial structure would exhibit neighbours throughout the domain. Fig. 1. View largeDownload slide An example RSN and its properties. The distance kernel is a Gaussian, $$f(d) = \exp(-d^2/2\sigma^2)/2\pi\sigma^2$$ with $$\sigma = 0.03$$. The imposed distribution of expected degrees is $$P(2)=P(15)=0.5$$. The density is $$\rho=10000$$. One node and its neighbours are highlighted. A random network without spatial structure would exhibit neighbours throughout the domain. A practical network model needs to be computationally feasible to generate. Appendix A describes an efficient generation algorithm for RSN networks, based on an $${\mathcal{O}}(N)$$ algorithm for Chung–Lu networks [44]. If the density is high enough that $$p_{uv}<1$$ for all node pairs, some important properties emerge. Consider a node $$u$$ at position $$\mathbf{x}$$ with a given $$\kappa_u$$, and consider network realizations containing $$u$$ but with different values of $$\rho$$. As $$\rho$$ increases, the distribution of the number of neighbours of $$u$$ in the network and the distribution of their locations does not change1. However, the probability that any two neighbours of $$u$$ are themselves neighbours of one another scales like $$1/\rho$$. Thus we anticipate that these networks will have similar spatial structure, but their clustering becomes negligible in the large $$\rho$$ limit. So the networks are locally tree-like at large $$\rho$$. This will allow our later analytic approaches to work. 2.1.1 Related network models We briefly review random network models that are related to RSNs. RSNs are a subclass of inhomogeneous random graphs [42, 45] and generalize many widely-studied models. They incorporate both degree heterogeneity and spatial structure in a way that, crucially, remains analytically tractable. They are a natural spatial extension of the well-known Erdős–Rényi and Chung–Lu random graphs. In a Chung–Lu network, an edge between $$u$$ and $$v$$ exists with probability $$\kappa_u\kappa_v/(N{\left\langle {\kappa} \right \rangle})$$ [23]. We recover this by setting $$f(d)=1/|V|$$ in Eq. (1). If we further set $$\kappa_u=\kappa$$ to be fixed for all $$u$$, we arrive at the “Erdős–Rényi” graph model, introduced by Gilbert [46]. RSNs are related to some other network models with spatial structure. There are a number of models which place nodes into a lattice and then assign edges by a rule that is at least in part determined by the distance between nodes. These include small-world networks [47, 48], navigable small-world networks [49], Waxman graphs [50], long-range percolation [51–53], scale-free percolation [54] (in particular [55] which uses a model very similar to our own) and spatially-constrained Erdős–Rényi graphs [56]. Most of the lattice-based models can be generated from the RSN model if we place the nodes on lattice points rather than uniformly at random. The widely-used random geometric graphs arise as a subclass of RSNs if we take uniform $$\kappa$$ and appropriately choose $$f$$. Furthermore, the class of spatially embedded random networks [57, 58], arises from a constant $$\kappa$$ and a power-law for $$f$$. The exponential random graph model does not have a clear connection to our model. It can handle degree heterogeneity [59] and some spatial structure [60], but typically the actual networks considered are small because they are expensive to generate. Many of these spatial networks are used to understand patterns emerging in the brain [18] or wireless sensor networks [61]. These models often have difficult-to-control correlations because of significant clustering. The microscopic structure makes analytic investigation of dynamic processes difficult. Due to this last weakness, almost all studies of spreading processes in networks with spatial structure are limited to simulation [62]. When equation-based analysis is attempted, it lacks quantitative agreement with simulations [56]. A major advantage of RSNs compared with other networks with spatial structure is their suitability to analytic results in the high-density limit. Although we demonstrate this for disease spread, many analytic techniques used to study other spreading processes in random networks will also apply to RSNs because of their locally tree-like structure in the high-density limit. We note however, that in some cases the high-density limit may not be achieved or significant clustering may exist, so the analytic model will not be accurate. In such cases, it may be possible to carry out an asymptotic expansion to identify corrections due to the lower density. We do not perform a detailed investigation of this. Throughout we assume that the nodes are embedded in a known geometry and a random network is generated for which the edges are assigned based on their distance and expected degrees. A closely related approach is given by latent space network models [63, 64]. In these models a network is given as the input, and the researcher attempts to find what underlying embedding geometry would most likely lead to such a network assuming a network generation algorithm similar to ours. There has been a significant amount of research based on this context, and many of these concepts also arise in studies on self-similar networks in hidden metric spaces [65]. Although much of the underlying mathematics in the network models is the same, our goal in using RSNs is to understand how processes would spread in a random network embedded in a given geometry, rather than to identify what geometry led to the creation of a given network. 2.2 Analytic governing equations for SIR disease on RSNs Our goal is to derive differential equations that accurately describe SIR disease dynamics on RSNs. We have noted that in the high-density limit RSNs have the crucial property of being locally tree-like, which allows us to apply and extend analytic tools that have been applied to other locally tree-like non-spatial networks like Chung–Lu or configuration model networks [35]. We derive analytic equations governing the deterministic dynamics (i.e., assuming that the number of nodes is sufficiently large for stochastic effects to become negligibly small). Our approach is of broad potential relevance, since we expect that analytic models on non-spatial graphs for processes other than SIR can be adapted to RSNs in a similar manner. In many random network classes, it is possible to develop powerful mathematical approaches by using their locally tree-like property: as $$N \to \infty$$ the clustering scales like $$1/N$$. This also occurs for our random spatial graphs: as $$\rho \to \infty$$, the clustering scales like $$1/\rho$$. We use this to develop analytic equations for SIR disease spread in the deterministic limit, following [35, 66, 67], but extended to the spatial setting. In an SIR epidemic, nodes begin susceptible and may be infected by infected neighbours (with rate $$\beta$$ per edge). Eventually infected nodes recover (with rate $$\gamma$$ per node). We assume that the population-scale dynamics are deterministic. That is, recognizing that the underlying process is stochastic, we assume that, as a function of $$\mathbf{x}$$, the proportion infected at some later time $$t$$ is uniquely determined from the initial proportion infected (as a function of $$\mathbf{x}$$). This assumption becomes reasonable in the limit of a large network, and is effectively the continuum assumption. Figure 2 illustrates that the spatio-temporal dynamics indeed become increasingly deterministic as the node density increases. Fig. 2. View largeDownload slide Disease spread in spatial networks in $$[0,1]\times[0.1]$$. The networks have $$P(3)=P(15)=1/2$$ and only short range connections [$$f(d)=1/0.01^2\pi$$ for $$d<0.01$$ and $$0$$ otherwise]. The value of $$\rho$$ varies, but has no effect on the bimodal degree distribution. Top: $$\rho=10^5$$, Middle: $$\rho=5\times 10^5$$, and Bottom: $$\rho=10^6$$. Disease is introduced to all nodes in a region of radius $$0.05$$ about the centre, with $$\beta=1/3$$ and $$\gamma=1$$. We plot the status of nodes at time $$t=12$$, with a detailed view of a subregion with the colours denoting a recovered node, an infected node or a susceptible node. As $$\rho$$ increases, the spread becomes more coherent and a deterministic limiting behaviour emerges. Fig. 2. View largeDownload slide Disease spread in spatial networks in $$[0,1]\times[0.1]$$. The networks have $$P(3)=P(15)=1/2$$ and only short range connections [$$f(d)=1/0.01^2\pi$$ for $$d<0.01$$ and $$0$$ otherwise]. The value of $$\rho$$ varies, but has no effect on the bimodal degree distribution. Top: $$\rho=10^5$$, Middle: $$\rho=5\times 10^5$$, and Bottom: $$\rho=10^6$$. Disease is introduced to all nodes in a region of radius $$0.05$$ about the centre, with $$\beta=1/3$$ and $$\gamma=1$$. We plot the status of nodes at time $$t=12$$, with a detailed view of a subregion with the colours denoting a recovered node, an infected node or a susceptible node. As $$\rho$$ increases, the spread becomes more coherent and a deterministic limiting behaviour emerges. We now derive analytic equations that faithfully describe the large $$\rho$$ SIR dynamics on RSNs. This will allow us to study and derive important properties of SIR disease propagation on networks with spatial structure. We start with the observation that the following four questions have identical answers in the large density limit, if indeed the population-scale dynamics are deterministic: As a function of location $$\mathbf{x}$$, what fraction of nodes in an infinitesimal region about $$\mathbf{x}$$ are susceptible, infected, or recovered at time $$t$$? What is the probability a random node is susceptible, infected, or recovered at time $$t$$, as a function of location $$\mathbf{x}$$? What is the probability a random node is susceptible, infected or recovered at time $$t$$, as a function of location $$\mathbf{x}$$, given the system’s state at time $$0$$? What is the probability a randomly chosen test node $$u$$ is susceptible, infected or recovered given the system’s state at time $$0$$ if we prevent $$u$$ from transmitting to its neighbours? The first two questions are clearly equivalent. The second and third are equivalent because we assume that the population-scale dynamics are deterministic. The third and fourth are equivalent because as long as $$u$$ is susceptible, the fact it does not transmit is irrelevant, and once $$u$$ is infected, its recovery time is not affected by any transmissions it causes. So the requirement that $$u$$ does not transmit does not influence the evolution of the state of $$u$$; it influences the states of neighbours of $$u$$, but only after it is too late for them to influence $$u$$. The process we analyse when preventing $$u$$ from transmitting (question 4) is different from the process in questions 1–3, but in question 4, we are not asking about the statuses of all nodes in the network. Rather we only ask about the status of node $$u$$, and for both processes, this is the same. Note that the equivalence of questions 3 and 4 holds specifically for SIR disease; it is crucial for our argument that once an individual begins to alter the status of its neighbours, its own future status can never be affected by those neighbours. We address question 4 using probabilistic tools. We define a test node $$u$$ which is randomly selected in the network at location $$\mathbf{x}$$ and prevented from transmitting to its neighbours. Under the assumption that $$\rho$$ is large, we can assume independence of $$u$$’s neighbours because clustering is negligible and because $$u$$ does not form a path of infection between its neighbours, precisely because we prevent it from transmitting. We seek to find the probability $$u$$ is susceptible, from which we deduce the probability it is infected or recovered. We pass to a continuum limit and define $$S(\mathbf{x},t)$$, $$I(\mathbf{x},t)$$ and $$R(\mathbf{x},t)$$ to be the probability that a test node $$u$$ placed at $$\mathbf{x}$$ would be susceptible, infected or recovered at time $$t$$. We similarly define $$\Phi_S(\mathbf{x},t)$$, $$\Phi_I(\mathbf{x},t)$$, $$\Phi_R(\mathbf{x},t)$$, and $$\Theta(\mathbf{x},t)$$ such that: $$\Phi_S$$ is the probability a random neighbour of the test node $$u$$ is susceptible at time $$t$$, $$\Phi_I$$ is the probability the neighbour is infected but has not transmitted to $$u$$, $$\Phi_R$$ is the probability the neighbour has recovered without transmitting to $$u$$ and $$\Theta=\Phi_S+\Phi_I+\Phi_R$$ is, thus, the probability that a random neighbour of node $$u$$ at position $$\mathbf{x}$$ has not transmitted infection to $$u$$ by time $$t$$. We assume infection is introduced to the population at time $$t=0$$ with $$S(\mathbf{x},0)$$ denoting the probability a node at $$\mathbf{x}$$ would be susceptible. We assume all other nodes are initially susceptible $$S(\mathbf{x},0) = 1- I(\mathbf{x},0)$$. Because $$k$$, the number of neighbours $$u$$ has, is Poisson distributed with mean $$\kappa_u$$ (the probability of a given $$k$$ is $$e^{-\kappa_u}\kappa_u^k/k!$$), the probability $$u$$ is susceptible at a later time $$t$$ is $$S(\mathbf{x},0) \sum_{k=0}^\infty e^{-\kappa_u}\kappa_u^k \Theta^k/k! = S(\mathbf{x},0) \exp[-\kappa_u(1-\Theta)]$$. Note that we have used here that the neighbours are independent (since we prevent $$u$$ from transmitting). We do not know the value of $$\kappa_u$$, so we must average over $$\kappa$$. The probability $$u$$ is susceptible is   \begin{align*} S &= S(\mathbf{x},0) \, \Psi(\Theta(\mathbf{x},t))\\ & = S(\mathbf{x},0) \int_0^\infty e^{-\kappa(1-\Theta(\mathbf{x},t))} P(\kappa) \, \mathrm{d}\kappa, \end{align*} with $$\Psi(\Theta(\mathbf{x},t)) \equiv \int_0^\infty e^{-\kappa(1-\Theta(\mathbf{x},t))} P(\kappa) \, \mathrm{d}\kappa$$. We augment this with $$I = 1-S-R$$ and $$\dot{R} = \gamma I$$. We must still find $$\Theta(\mathbf{x},t)$$. We turn to the flow diagram in Fig. 3. Once a neighbour $$v$$ of the test node $$u$$ becomes infected, whether or not $$v$$ transmits infection to $$u$$ and if so, at what time the transmission occurs, is independent of anything else that happens in the population. Thus we can immediately calculate that the flux from $$\Phi_I$$ to $$1-\Theta$$ is $$\beta\Phi_I$$, so $$\dot{\Theta} = -\beta \Phi_I$$. Similarly we find that the flux to $$\Phi_R$$ is $$\gamma \Phi_I$$. Since the flux to $$1-\Theta$$ and $$\Phi_R$$ are both proportional to $$\Phi_I$$ and both are $$0$$ at $$t=0$$, we can conclude that $$\Phi_R = \gamma(1-\Theta)/\beta$$. However, the rate at which $$v$$ becomes infected is more difficult. It will be easier to directly calculate $$\Phi_S$$ in terms of $$\Theta$$ and use $$\Phi_I = \Theta-\Phi_S-\Phi_R$$ to give an equation for $$\dot{\Theta}$$ in terms of $$\Theta$$. Fig. 3. View largeDownload slide A flow diagram leading to the governing equation for $$\Theta(\mathbf{x},t)$$. The compartments $$\Phi_S$$, $$\Phi_I$$ and $$\Phi_R$$ make up $$\Theta$$ and represent the probability that a random partner has not transmitted to $$u$$ and is susceptible, infected or recovered, respectively. Fig. 3. View largeDownload slide A flow diagram leading to the governing equation for $$\Theta(\mathbf{x},t)$$. The compartments $$\Phi_S$$, $$\Phi_I$$ and $$\Phi_R$$ make up $$\Theta$$ and represent the probability that a random partner has not transmitted to $$u$$ and is susceptible, infected or recovered, respectively. To find $$\Phi_S$$, we consider the possible neighbours of $$u$$ and calculate their probability of being susceptible. The probability a node $$v$$ at $$\hat{\mathbf{x}}$$ is a neighbour of $$u$$ is proportional to $$\kappa_v$$ and to $$f(|\hat{\mathbf{x}}-\mathbf{x}|)$$. In turn, the probability $$v$$ is susceptible is $$S(\hat{\mathbf{x}},0) \exp[-\kappa_v(1-\Theta(\hat{\mathbf{x}},t))]$$. So the probability a random neighbour is susceptible is given by   \begin{align*} \Phi_S &= \int_V S(\hat{\mathbf{x}},0)\int_\kappa \frac{\kappa P(\kappa)}{{\left\langle {\kappa} \right \rangle}} f(|\hat{\mathbf{x}}-\mathbf{x}|) e^{-\kappa(1-\Theta(\hat{\mathbf{x}},t))} \, \mathrm{d}\kappa \, \mathrm{d}\hat{\mathbf{x}}\\ &= \frac{\int_V S(\hat{\mathbf{x}},0) \Psi'(\Theta(\hat{\mathbf{x}},t)) f(|\hat{\mathbf{x}}-\mathbf{x}|) \, \mathrm{d}\hat{\mathbf{x}}}{{\left\langle {\kappa} \right \rangle}}. \end{align*} From this, $$\frac{\partial {}}{\partial {t}}{\Theta} = -\beta\Phi_I=-\beta(\Theta-\Phi_R-\Phi_S)$$ becomes   \begin{align*} \frac{\partial {}}{\partial {t}}{\Theta}(\mathbf{x},t) &= - \beta \Theta(\mathbf{x},t) + \gamma(1-\Theta(\mathbf{x},t))\\ &\quad + \beta \frac{\int_{V}S(\hat{\mathbf{x}},0) \Psi'(\Theta(\hat{\mathbf{x}},t)) f(|\hat{\mathbf{x}}-\mathbf{x}|) \, \mathrm{d}\hat{\mathbf{x}}}{{\left\langle {\kappa} \right \rangle}}. \end{align*} This derivation yields the following set of coupled equations for SIR propagation on RSNs:   \begin{align} \frac{\partial {}}{\partial {t}}{\Theta}(\mathbf{x},t) &= - \beta \Theta(\mathbf{x},t) + \gamma(1-\Theta(\mathbf{x},t)) \nonumber \label{eq:dtheta}\\ &\quad+ \beta \frac{\int_{V}S(\hat{\mathbf{x}},0)\Psi'(\Theta(\hat{\mathbf{x}},t)) f(|\hat{\mathbf{x}}-\mathbf{x}|) \, \mathrm{d}\hat{\mathbf{x}}}{{\left\langle {\kappa} \right \rangle}}\\ \end{align} (2a)  \begin{align} S(\mathbf{x},t) &= S(\mathbf{x},0)\Psi(\Theta(\mathbf{x},t)),\\ \end{align} (2b)  \begin{align} \frac{\partial {}}{\partial {t}}{R}(\mathbf{x},t) &= \gamma (1-S(\mathbf{x},t)-R(\mathbf{x},t))\\ \end{align} (2c)  \begin{align} I(\mathbf{x},t) &=1-S(\mathbf{x},t)-R(\mathbf{x},t) \label{eq:I} \end{align} (2d) with $$\Psi(\Theta) = \int_0^\infty e^{-\kappa(1-\Theta)} P(\kappa) \, \mathrm{d}\kappa$$. Our initial conditions are $$\Theta(\mathbf{x},0)=1$$ and $$R(\mathbf{x},0)=0$$, with $$S(\mathbf{x},0) = 1- I(\mathbf{x},0)$$. The equation for $$\Theta$$ is a nonlocal evolution equation with an integral over space and a derivative with respect to time. The nonlocal effects are captured through the convolution integral. To emphasize that the nonlocal interactions are captured by an integral, we refer to this as a partial integro-differential equation (PIDE): Unlike [35] we obtain an evolution equation for $${\Theta}(\mathbf{x},t)$$ that captures the RSN’s spatial structure in the integral term. We expect the equation for $$\Theta$$ to be similar to the Fisher–KPP equation for a spreading population [68, 69] $$\dot{u} = ru(1-u/K) + D u_{xx}$$, with the spatial integral playing the role of the nonlinear and diffusion terms. The spatial integral captures the network’s structure by including nonlocal interactions through $$f$$ and the degree distribution through $$\Psi$$. In Appendix B, we derive the Fisher–KPP equation as an approximation of (2a). Other nonlocal versions of the Fisher–KPP equation have been studied [70, 71], including some for disease spread [9, 39–41]. These are based on a mass-action mixing assumption and involve a convolution of a distance kernel with $$u$$ directly rather than a nonlinear function of $$u$$. Since stochastic simulations on networks are often hard to analyse and understand, this explicit analytic equation provides a powerful tool to study SIR dynamics on random networks with versatile spatial structure, fully taking into account the distribution of expected degrees and the distance kernel. In particular, the stochastic simulations to be discussed in Section 3 reveal the existence of travelling waves. In the next subsection, we use the PIDE to analytically derive properties of travelling waves in the analytic model, which, as we demonstrate in Section 3, accurately describe the travelling waves observed in the stochastic simulations. 2.3 Analytic properties of travelling waves on RSNs The PIDE formulation provides a powerful technical tool to derive precise quantitative insight in the spread of SIR disease on RSNs. As a major application, we can derive properties of travelling waves allowed by the PIDE. In particular, we derive an explicit expression for the wave speed of a 1D travelling wave, and identify conditions on the spatial decay of the kernel for the travelling wave to exist. This wave has a ‘pulled front’: That is, its speed is set by the nodes in the leading edge. We derive the wave speed assuming that the domain $$V$$ is the real line. We can write $$\Theta(x,t) = \theta(\xi(x,t))$$ where $$\xi(x,t) = x-ct$$ and $$c$$ is the wave speed. We will assume it is travelling from left to right. The wave travels into a region where $$S(x,t)\approx 1$$, $$I(x,t) \ll 1$$ and $$R(x,t) \ll 1$$. Very little transmission has occurred, so ahead of the travelling wave $$\Theta \approx 1$$. We assume $$\xi$$ is in the leading edge and write $$\Theta(x,t) = \theta(\xi(x,t))=1-\epsilon h(\xi(x,t))$$. In this leading edge of the wave $$\epsilon h(\xi) \ll 1$$, while in the bulk of the wave, $$\epsilon h(\xi)$$ may be comparable to $$1$$. We focus on the behaviour in the leading edge of the wave and transform equation (2a) for $$\frac{\partial {}}{\partial {t}}\Theta$$ into an equation for $$h$$ by expanding $$\Psi'(\theta(\xi))$$ as a Taylor Series about $$\theta=1$$:   \[ \Psi'(\theta) = \Psi'(1) - \epsilon h(\xi) \Psi''(1) + \frac{\epsilon^2 h(\xi)^2}{2}\Psi'''(1) + \cdots \] Note that $$\Psi^{(n)}(1) = {\left\langle {\kappa^n} \right \rangle}$$. Substituting this into the integral (and taking $$S(x,0)=1$$ ahead of the wave), we arrive at   \begin{align*} \epsilon c h'(\xi) &= -\beta(1-\epsilon h(\xi)) + \epsilon \gamma h(\xi) \\ &\quad +\beta \frac{\int_{-\infty}^\infty [\Psi'(1) - \epsilon h(\hat{\xi}) \Psi''(1)]\,f(|\hat{\xi}-\xi|) \, \mathrm{d}\hat{\xi}}{{\left\langle {\kappa} \right \rangle}}\\ & \quad + \beta \frac{\int_{-\infty}^\infty{\mathcal{O}}(\epsilon^2h(\hat{\xi})^2)\Psi'''(1)\,f(|\hat{\xi}-\xi|) \mathrm{d}\hat{\xi}}{{\left\langle {\kappa} \right \rangle}} \end{align*} Because $$\Psi'(1)={\left\langle {\kappa} \right \rangle}$$, the $${\mathcal{O}}(1)$$ terms cancel. We neglect the $${\mathcal{O}}(\epsilon^2 h(\hat{\xi})^2)$$ terms by arguing that if $$|\hat{\xi}-\xi|$$ is not large then $$\epsilon^2h(\hat{\xi})^2$$ is small (at the leading edge), and if $$|\hat{\xi}-\xi|$$ is large then $$f(|\hat{\xi}-\xi|)$$ is small (away from the leading edge). This leaves the linear homogeneous equation for $$h$$  \[ ch'(\xi) = (\beta+\gamma)h(\xi) - \beta \frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}\int_{-\infty}^\infty h(\hat{\xi})f(|\hat{\xi}-\xi|) \, \mathrm{d}\hat{\xi} \] We anticipate $$h(\xi) \approx e^{-\alpha \, \xi}$$ for some $$\alpha$$, yielding:   \[ -c\alpha = (\beta+\gamma) - \beta \frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}\int_{-\infty}^\infty e^{-\alpha (\hat{\xi}-\xi)} f(|\hat{\xi}-\xi|) \, \mathrm{d}\hat{\xi} \] Setting $$u = \hat{\xi}-\xi$$ the integral becomes $$\int_{-\infty}^\infty e^{-\alpha u} f(|u|)\, \mathrm{d}u$$. This is the bilateral Laplace transform of $$f(|x|)$$, which we denote $${\mathcal{L}}{}[\,f](\alpha)$$. Performing some further algebra yields   \begin{equation} \frac{c}{\beta+\gamma} = -\frac{1}{\alpha} + {\mathcal{R}_0} \frac{{\mathcal{L}}{}[\,f](\alpha)}{\alpha} \label{eq:alpha_c} \end{equation} (3) where $${\mathcal{R}_0} = \beta{\left\langle {\kappa^2} \right \rangle}/(\beta+\gamma){\left\langle {\kappa} \right \rangle}$$ (this is the basic reproductive number of the SIR disease on the network, which is the typical number of infections caused by an infected individual early in an epidemic [31]). There are infinitely many solution pairs $$\alpha,c$$. Following results for the Fisher–KPP equation we expect that the observed solution has the minimum wave speed. Setting $$\frac{\partial {c}}{\partial {\alpha}}=0$$ and performing some algebra shows that this occurs when   \begin{equation} \alpha{\mathcal{L}}{}[xf(|x|)](\alpha) + {\mathcal{L}}{}[\,f](\alpha) = \frac{1}{{\mathcal{R}_0}} \label{eqn:find_alpha} \end{equation} (4) We can solve this equation to find $$\alpha$$. Finally, substituting (3) into (4) gives   \begin{equation} c = -\beta\frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}} {\mathcal{L}}{}[xf(|x|)](\alpha) = -(\beta+\gamma) {\mathcal{R}_0} {\mathcal{L}}{}[xf(|x|)](\alpha) \label{eqn:find_c} \end{equation} (5) Note that the the wave speed $$c$$ and the spatial decay rate $$\alpha$$ depend on the distribution of expected degrees only through $${\left\langle {\kappa} \right \rangle}$$ and $${\left\langle {\kappa^2} \right \rangle}$$, specifically through the ratio $${\left\langle {\kappa^2} \right \rangle}/{\left\langle {\kappa} \right \rangle}$$. Figure 4 investigates the dependence of the speed on the degree distribution. Specifically, for different $$f$$ proportional to $$e^{-|x|^{1+\epsilon}}$$, we investigate how the wavespeed depends on $${\mathcal{R}_0} = \beta {\left\langle {\kappa^2} \right \rangle}/{\left\langle {\kappa} \right \rangle}$$. We normalize time so that $$\beta+\gamma=1$$. As $${\mathcal{R}_0}$$ grows, eventually $$\alpha$$ changes very slowly, and $$c$$ is approximately proportional to $${\mathcal{R}_0}$$. So increasing $${\mathcal{R}_0}$$ by modifying the distribution of $$\kappa$$ increases the wavespeed, and the speed tends to be proportional to $${\left\langle {\kappa^2} \right \rangle}/{\left\langle {\kappa} \right \rangle}$$. Fig. 4. View largeDownload slide Dependence of predicted wave properties on degree distribution. The predicted values of $$\alpha$$ and $$c$$ as a function of $${\mathcal{R}_0}$$ for $$f(x) \propto e^{-|x|^{1+\epsilon}}$$. We measure how the wavespeed depends on $${\mathcal{R}_0} = \frac{\beta}{\beta+\gamma}\frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}$$ with time normalized so that $$\beta+\gamma=1$$. Fig. 4. View largeDownload slide Dependence of predicted wave properties on degree distribution. The predicted values of $$\alpha$$ and $$c$$ as a function of $${\mathcal{R}_0}$$ for $$f(x) \propto e^{-|x|^{1+\epsilon}}$$. We measure how the wavespeed depends on $${\mathcal{R}_0} = \frac{\beta}{\beta+\gamma}\frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}$$ with time normalized so that $$\beta+\gamma=1$$. Note that if we chose another normalization of time, such as setting $$\gamma=1$$ our results suggest that increasing $${\mathcal{R}_0}$$ by increasing $$\beta$$ would mean that the wavespeed would scale like $$(1+\beta)\beta$$ for larger $$\beta$$. Note also that at $${\mathcal{R}_0}=1$$, both $$\alpha$$ and $$c$$ are $$0$$. If $$f$$ does not decay at least exponentially fast, then the Laplace transforms above do not exist. Thus we infer that if $$f$$ does not decay at least exponentially fast, then there is no coherent travelling wave solution in the $$\rho\to\infty$$ limit, consistent with [41, 72]. In practice, for a finite population if the long tail is observed by the transmissions, we see hop-and-spread dynamics, while if the long tail is not sampled by the transmission we may still see travelling wave behaviours in a given realization. This result for the propagation of SIR disease on random spatial networks described by Eq. (1) precisely reproduces the seminal result of Mollison [41] that for a simple 1D model of disease propagation in continuous space, travelling wave solutions can exist only for kernels that decay exponentially or faster. This is no surprise, because although the model used by Mollison (and, in a similar fashion, in [9]), was formulated in a simplified setting involving a convolution of a distance kernel with $$u$$ directly rather than a nonlinear function of $$u$$, it is still similar to our new PIDE. It is a powerful result, though, that we can analytically demonstrate this same behaviour, originally derived for simplified continuous space models, on the discrete networks generated by the RSN formula Eq. (1), in the limit of large node density. 2.4 Epidemic final size A further prediction for stochastic SIR spread on RSNs that can be computed analytically from the PIDE is the final size of an epidemic in a large population. As $$t \to \infty$$ the system approaches an equilibrium. By setting $$\frac{\partial {}}{\partial {t}}\Theta = 0$$, we have   \[ \Theta(\mathbf{x},\infty) = \frac{\gamma}{\beta+\gamma} + \frac{\beta}{\beta+\gamma} \frac{S(\hat{\mathbf{x}},0)\int_V \Psi'(\Theta(\hat{\mathbf{x}},\infty)) f(|\hat{\mathbf{x}}-\mathbf{x}|) \, \mathrm{d}\hat{\mathbf{x}}}{{\left\langle {\kappa} \right \rangle}} \] If $$\mathbf{x}$$ is far from the point of introduction, or the introduction is so small that we can approximate $$S(\mathbf{x},0)$$ as $$1$$ everywhere, we can treat $$\Theta(\mathbf{x},\infty)$$ as spatially homogeneous. Then   \begin{align} \Theta &= \frac{\gamma}{\beta+\gamma} + \frac{\beta}{\beta+\gamma} \frac{\Psi'(\Theta)}{{\left\langle {\kappa} \right \rangle}} \, , \label{eqn:theta_fs}\\ \end{align} (6)  \begin{align} S &= \Psi(\Theta) \, . \end{align} (7) $$\Theta=1$$, $$S=1$$ is always a solution, but if there is a solution $$\Theta$$ between $$0$$ and $$1$$, it is the appropriate one for an epidemic. It exists precisely when $${\mathcal{R}_0}>1$$. This is the same relation as for a random Chung–Lu network without spatial structure [35]. Interestingly, this means the epidemic final size in RSNs is independent of the distance kernel and depends only on disease properties and the degree distribution. 3. Numerical results We now apply the analytic models from Section 2 to analyse and understand stochastic simulations of SIR disease spread on RSNs. We first consider a simple Gaussian spatial kernel and demonstrate that the nonlocal evolution equations quantitatively predict stochastic simulations. The theory of Section 2 predicts the occurrence of travelling waves for the Gaussian kernel, and we do indeed observe travelling waves in the stochastic simulations on RSNs that are closely matched by numerical solutions of travelling waves in the PIDE. We also verify the analytic properties that were derived for the travelling waves in Section 2. We then consider a specific subclass of RSNs on which we investigate the interplay between short-range and long-range spatial connections in the SIR dynamics. 3.1 travelling waves and correspondence between stochastic SIR simulation and PIDE model We first compare numerical simulations of PIDE System (2) with stochastic SIR simulations on RSNs with a simple Gaussian kernel to demonstrate that these equations closely describe the stochastic dynamics in the limit of large $$N$$. In the process, we will identify travelling waves and verify the analytic properties we derived for them. We consider stochastic simulations on one-dimensional (1D) and two-dimensional (2D) RSNs generated using an algorithm based on the linear-time algorithm of [73] for Chung–Lu networks (see Appendix A). We simulate epidemic spread starting from a localized initial condition using a stochastic simulation algorithm from [32]. Figure 5 shows travelling waves revealed by the stochastic simulations in 1D and 2D. These retain their shape as they propagate. Figure 5 also shows numerical PIDE solutions of System (2) for the 1D and 2D network problems, demonstrating good agreement between PIDE solution and stochastic simulations on the 1D and 2D RSNs with Gaussian spatial kernel. Fig. 5. View largeDownload slide Comparison of stochastic simulation of SIR disease dynamics with numerical solution of analytic PIDE, on 1D and 2D RSNs. In both panels, results of stochastic simulation on networks with $$N=10^6$$ nodes are given by black dots. (Left) 1D PIDE solution is given by solid blue line (spatial resolution $$\Delta x = 10^{-4}$$). (Right) 2D PIDE solution is given by mesh surface (spatial resolution $$\Delta x = 1/512$$). Networks in both panels are generated using kernel $$f(r) = \phi_{0,0.01}(r)$$, where $$\phi_{\mu,\sigma}(x)$$ is the probability density function for the normal random variable with mean $$\mu$$ and standard deviation $$\sigma$$. All nodes have expected degree $$\kappa=20$$. Initial conditions: (left) 10% of nodes in the interval $$[0,0.01)$$ are initially infected, and (right) 10% of nodes in the square $$[0,1/32)\times[0,1/32)$$ are initially infected. Fig. 5. View largeDownload slide Comparison of stochastic simulation of SIR disease dynamics with numerical solution of analytic PIDE, on 1D and 2D RSNs. In both panels, results of stochastic simulation on networks with $$N=10^6$$ nodes are given by black dots. (Left) 1D PIDE solution is given by solid blue line (spatial resolution $$\Delta x = 10^{-4}$$). (Right) 2D PIDE solution is given by mesh surface (spatial resolution $$\Delta x = 1/512$$). Networks in both panels are generated using kernel $$f(r) = \phi_{0,0.01}(r)$$, where $$\phi_{\mu,\sigma}(x)$$ is the probability density function for the normal random variable with mean $$\mu$$ and standard deviation $$\sigma$$. All nodes have expected degree $$\kappa=20$$. Initial conditions: (left) 10% of nodes in the interval $$[0,0.01)$$ are initially infected, and (right) 10% of nodes in the square $$[0,1/32)\times[0,1/32)$$ are initially infected. Figure 6 compares the analytically predicted wave speed $$c$$ with wave speeds obtained from stochastic simulations. The small black circles of Fig. 6 give the wave speed observed in stochastic simulations. They converge to the analytic prediction (cyan dashed), but the convergence is very slow. Due to finite density, the exponentially decaying front is eventually truncated at its leading edge as shown in Fig. 7. This truncation slows the wave because the foremost infected nodes play a large role in setting the wave speed. This has been analysed for similar fronts in other stochastic systems, for which the leading order error in the wave speed decays proportionally to $$1/[\ln \rho]^2$$ [74]. More nodes are needed to improve the fit, as seen in Fig. 6. Figure 7 confirms that the stochastic simulation front and the numerically calculated PIDE front are nearly exponential with slope close to the predicted $$-\alpha$$. For the stochastic simulation (green dots), local densities in the exponentially decaying profile smaller than about 10$$^{-4}$$ cannot be represented because there are insufficient points in the simulation (the local densities effectively drop down to zero to the right of $$x\sim 0.5$$, and these zero values end up outside the range of the vertical logarithmic axis of the figure). By having larger $$\rho$$, more of the leading edge is observed resulting in wave speeds closer to the analytic prediction (Fig. 6). Figure 7 also shows that, in the numerical PIDE simulations, the exponential profile is truncated when the density of infected individuals approaches machine accuracy. Fig. 6. View largeDownload slide Comparison of wave speed for SIR disease dynamics observed in stochastic simulations (small black circles) and analytic prediction (dashed cyan line). For the stochastic simulations, we generate $$n_{rep} = 25$$ 1D spatial networks for each node density $$\rho$$ using a distance kernel $$f(|x|) = \phi_{0,0.01}(|x|)$$, where $$\phi_{\mu,\sigma}(x)$$ is the probability density function for the normal random variable with mean $$\mu$$ and standard deviation $$\sigma$$. All nodes have expected degree $$\kappa=10$$. For each network we realize one SIR simulation with disease parameters $$\beta=1$$ and $$\gamma=3$$. The black circles show the average wave speed resulting from 25 network realizations (vertical bars indicate 95% confidence intervals). The dotted black curve represents the expected convergence behaviour of $$c^*- K/[\ln \rho]^2$$ [74] to the analytically predicted wave speed $$c^*$$, where the constant $$K$$ is obtained by fitting the curve to the rightmost black circle, which has the smallest error bar since it corresponds to the largest node density $$\rho$$. Fig. 6. View largeDownload slide Comparison of wave speed for SIR disease dynamics observed in stochastic simulations (small black circles) and analytic prediction (dashed cyan line). For the stochastic simulations, we generate $$n_{rep} = 25$$ 1D spatial networks for each node density $$\rho$$ using a distance kernel $$f(|x|) = \phi_{0,0.01}(|x|)$$, where $$\phi_{\mu,\sigma}(x)$$ is the probability density function for the normal random variable with mean $$\mu$$ and standard deviation $$\sigma$$. All nodes have expected degree $$\kappa=10$$. For each network we realize one SIR simulation with disease parameters $$\beta=1$$ and $$\gamma=3$$. The black circles show the average wave speed resulting from 25 network realizations (vertical bars indicate 95% confidence intervals). The dotted black curve represents the expected convergence behaviour of $$c^*- K/[\ln \rho]^2$$ [74] to the analytically predicted wave speed $$c^*$$, where the constant $$K$$ is obtained by fitting the curve to the rightmost black circle, which has the smallest error bar since it corresponds to the largest node density $$\rho$$. Fig. 7. View largeDownload slide A comparison of stochastic simulation and numerical PIDE solution on a logarithmic scale, using the 1D solutions from Fig. 5. Both stochastic and numerical solutions experience leading edge truncation. For the stochastic simulation it is due to the finite number of nodes, while for the numerical solution it is due to numerical precision ($$\sim 10^{-16}$$). The slope of the leading edge is close to the asymptotic prediction of $$-\alpha$$. Fig. 7. View largeDownload slide A comparison of stochastic simulation and numerical PIDE solution on a logarithmic scale, using the 1D solutions from Fig. 5. Both stochastic and numerical solutions experience leading edge truncation. For the stochastic simulation it is due to the finite number of nodes, while for the numerical solution it is due to numerical precision ($$\sim 10^{-16}$$). The slope of the leading edge is close to the asymptotic prediction of $$-\alpha$$. 3.2 A subclass of RSNs with long-range connections We now consider a specific subclass of RSNs with short-range and long-range spatial connections on which we will investigate the effect long-range connections have on SIR dynamics. The networks in this subclass have an average degree of $$k$$, with edges distributed such that nodes within some distance $$r_0$$ are very likely to be connected, and nodes of greater distance are very unlikely to be connected. These are modelled after the Watts–Strogatz small-world networks [47] with the number of local edges decreasing as the long-range connections increase so that the average degree remains fixed. Specifically, we set $$\kappa = k$$ for all nodes, choose $$r_0$$, and define a reference node density $$N_0:=k/(\pi r_0^2)$$. We place $$N \ge N_0$$ nodes uniformly at random into the unit torus $$[0,1]\times[0,1]$$ with periodic boundaries and $$|V|=1$$. We assume $$r_0<1/2$$ so that disks of radius $$r_0$$ easily fit within the torus. In particular this forces $$\pi r_0^2<1$$. We choose a small number $$\epsilon$$ and define the spatial kernel as   \begin{equation} f(d_{uv}) = \begin{cases} \frac{N_0}{k}\left[1-\epsilon\frac{1-\pi r_0^2}{\pi r_0^2}\right]& d_{uv} < r_0\\ \frac{N_0}{k}\epsilon & d_{uv} \geq r_0 \end{cases}. \label{eq:kernel} \end{equation} (8) Note that $$\epsilon$$ cannot be greater than $$\pi r_0^2/(1-\pi r_0^2)$$. Equation (1) specializes to   \begin{equation} p_{uv} = \min\left(k \frac{f(d_{uv})}{N}, 1\right) = \frac{k f(d_{uv})}{N}, \label{eq:pspec} \end{equation} (9) since, for small $$\epsilon$$, $$f(d_{uv}) \le N_0/k$$ so $$kf(d_{uv})/N \le 1$$, assuming $$\pi r_0^2<1$$. Then   \begin{equation} p_{uv} = \begin{cases} \frac{N_0}{N} \left(1- \epsilon\frac{1-\pi r_0^2}{\pi r_0^2}\right) & d_{uv} <r_0\\ \frac{N_0}{N} \epsilon & d_{uv}\geq r_0 \end{cases}. \label{eq:puv} \end{equation} (10) The expected number of nodes within distance $$r_0$$ is $$\pi r_0^2 N$$ and beyond distance $$r_0$$ is $$(1-\pi r_0^2)N$$. The total number of expected neighbours in the graph is $$\pi r_0^2 N \frac{N_0}{N} \left(1-\epsilon(1-\pi r_0^2)/(\pi r_0^2)\right)+(1-\pi r_0^2)N \frac{N_0}{N} \epsilon=k$$. When we place $$N=N_0$$ nodes and set $$\epsilon=0$$, all nodes within a distance $$r_0$$ are connected ($$p_{uv}=1$$) and no other nodes are joined. This is a random geometric graph [75]. When $$\epsilon=0$$ and $$N$$ is chosen greater than $$N_0$$, there are more than $$k$$ nodes within a distance $$r_0$$, but not all of these nodes are connected ($$p_{uv} < 1$$), such that the expected degree is still $$k$$. For any choice of $$N \geq N_0$$, as $$\epsilon$$ increases, the number of long-range connections grows proportionally to $$\epsilon$$ with a corresponding decrease in the short-range connections. In what follows we consider networks of this subclass with $$N>N_0$$ and varying $$\epsilon$$, i.e., a varying fraction of long-range connections. 3.3 An example of SIR disease propagation on RSNs with short-range and long-range connections: travelling waves and hop-and-spread dynamics We now perform stochastic simulations [32] of SIR disease on RSNs with spatial kernel (8) and $$k=20$$, for three values of $$\epsilon$$. We choose $$N_0= 10^5$$ such that $$r_0=\sqrt{k/(\pi N_0)} = 8.0\times10^{-3}$$, and place $$N=10^6$$ random nodes. Nodes transmit infection with rate $$\beta=1$$ and recover with rate $$\gamma=3$$. The simulated behaviour shown in Fig. 8 reveals: When there are few long-range connections ($$\epsilon \to 0$$), the disease spreads in a wave-like, coherent spatial pattern outwards from the source. With an increased number of long-range connections, we see hop-and-spread dynamics. The spread is locally wave-like until a long-range transmission occurs. It then generates a new focal region. Once the number of infections is large enough that the disease experiences its first long-range transmission, it usually starts causing many long-range transmissions. As the number of long-range transmissions increases, the disease spreads uniformly throughout the population. Sufficient long-range connections mean the epidemic behaves like it would in a globally-mixed population. Fig. 8. View largeDownload slide Stochastic simulation of SIR disease dynamics on a 2D RSN with spatial kernel (8) with $$N=10^6$$ nodes. The fraction of infected nodes, $$I$$, is shown, with a colour scale that ranges from blue ($$I=0$$) to yellow ($$I=0.5$$). The density of infected nodes at time $$t=10$$ seconds is shown for networks with (left) highly local spatial structure ($$\epsilon=10^{-10}$$), (middle) intermediate spatial structure ($$\epsilon=10^{-8.25}$$) and (right) no local spatial structure ($$\epsilon=\pi r_0^2$$, i.e., $$f(d_{uv})=1$$, uniform spatial kernel). Initial condition: the 1000 nodes closest to (0.5,0.5) are initially infected. The waves in the first two panels behave as travelling waves, as can be seen in the supplementary movies 1–3 (one for each panel). Fig. 8. View largeDownload slide Stochastic simulation of SIR disease dynamics on a 2D RSN with spatial kernel (8) with $$N=10^6$$ nodes. The fraction of infected nodes, $$I$$, is shown, with a colour scale that ranges from blue ($$I=0$$) to yellow ($$I=0.5$$). The density of infected nodes at time $$t=10$$ seconds is shown for networks with (left) highly local spatial structure ($$\epsilon=10^{-10}$$), (middle) intermediate spatial structure ($$\epsilon=10^{-8.25}$$) and (right) no local spatial structure ($$\epsilon=\pi r_0^2$$, i.e., $$f(d_{uv})=1$$, uniform spatial kernel). Initial condition: the 1000 nodes closest to (0.5,0.5) are initially infected. The waves in the first two panels behave as travelling waves, as can be seen in the supplementary movies 1–3 (one for each panel). The coherent spatial spread and globally-mixed regimes behave largely deterministically, but the intermediate hop-and-spread regime is dominated by stochastic effects. In a large enough domain (or equivalently, for small enough $$r_0$$), an epidemic spreading with any $$\epsilon>0$$ will eventually exhibit hop-and-spread dynamics. How long the initial wave travels before the hop-and-spread dynamics begin depends on the time until the disease crosses a long-range connection. We can gain insight by initially ignoring the long-range connections. Then the disease spreads in an annulus. The width and rate of spread of the annulus is roughly constant as the radius increases. Thus the number of infected nodes increases linearly. From the frequency of long-range edges and the transmission probability per edge, we can estimate the expected number of infections until a long-range transmission occurs. Once one long-range transmission occurs, the number infected is likely high enough that many other long-range transmissions begin to occur. So the typical dynamics have an early phase in which local transmissions dominate, yielding a growing localized annulus of infection. Then long-range transmissions followed by new foci occur, and this dominates the later dynamics. 3.4 Factors affecting accuracy of analytic approximation The analytic equations are effectively a continuum approximation to the discrete dynamics of the process. To be a valid approximation, it needs to be sensible to talk about the density of nodes and edges near a point in space. That is, we need the transition from Question 1 to Question 2 in section 2.2 to be appropriate. Thus we need a separation of scales: we need a scale which is small compared with the dynamics of the spatial spread, but large enough that the local expected behaviour gives a very accurate prediction of the number of nodes within the region having a given status. The specific conditions under which the analytic equations become appropriate are both disease and network-dependent. Let us consider a region $$X$$ which is small compared with the spatial scales over which the disease spread is predicted to vary. The predictions will tend to be inaccurate when they predict a very small number of infected individuals in $$X$$. We can classify these into four situations If the local reproductive number is less than one, the inaccuracies are not important because the disease is dying off in both the prediction and the realization. If $$X$$ is slightly ahead of an advancing front, the error will tend to cause a reduction in the front’s speed, but otherwise the predictions perform well. If $$X$$ is far ahead of any front, then the impact may be significant: an introduced single case where the expected number of cases is much less than $$1$$ can lead to a new focus of spread, as seen in Fig. 8. If the density of nodes is low enough that $$X$$ does not have many nodes or clustering is significant, then the arguments leading to the governing equations break down. The first two issues lead to relatively small errors, while the final two can lead to significant errors. 4. Discussion and conclusion In this article, we have derived new analytic models of PIDE type for the spread of SIR disease on random networks with spatial structure. We derive these models for a class of RSNs that intrinsically incorporates spatial structure in a way that crucially remains analytically tractable. These models are important to understand and quantify the dynamics of SIR disease spread on networks with spatial structure. We observe nonlinear travelling waves on RSNs and are able to precisely characterize these waves analytically. This is the first quantitatively accurate analytic description of nonlinear travelling waves on random graphs with spatial structure. Our SIR disease model on spatial graphs extends analytic approaches that were previously only available to describe dynamics at the level of populations or on networks without spatial structure, to the real-life case of networks characterized by important spatial structure. This provides a theoretical modelling framework, for example, for recent observations of how epidemics like Ebola evolve in modern connected societies, with long-range connections seeding new focal points from which the epidemic spreads locally (see [11], in particular Video 1). Figure 9 shows the observed hops seen from virus genome sequencing for the 2013-2016 Ebola outbreak [11], which can be built into RSN models. Fig. 9. View largeDownload slide Observed hop distances inferred from viral genome sequencing for the 2013–2016 Ebola outbreak in West-Africa. Data courtesy of Rambaut [11]. The estimated probability density is obtained from raw binned data by Kernel Density Estimation using a Gaussian kernel function. The dotted lines denote the 50th and 95th percentiles. Fig. 9. View largeDownload slide Observed hop distances inferred from viral genome sequencing for the 2013–2016 Ebola outbreak in West-Africa. Data courtesy of Rambaut [11]. The estimated probability density is obtained from raw binned data by Kernel Density Estimation using a Gaussian kernel function. The dotted lines denote the 50th and 95th percentiles. There are many possibilities for further applications in disease modelling. For example, the SIR travelling wave models can be applied to realistic simulations of historic epidemics such as the plague in medieval Europe [9], incorporating in a precise manner the geography and estimated historical population density maps. Further potential applications include models for a diversity of areas such as neuronal networks in the human brain, spread of animal and plant species in fragmented habitats, and propagation of civil unrest, or public health epidemics such as obesity and smoking, in human societies. The RSN class is also amenable to graph theoretic analysis. For example, it would be of great interest to investigate thresholds for the existence of a giant component for RSNs with various distance kernels. This was successfully done for non-spatial random graphs with given degree sequences by Reed et al. [36, 76]. The techniques in [77] show that this should be straightforward for RSNs. Also, the ODE systems resulting from applying the analytic edge-based compartmental modelling to SIR disease spread on networks without spatial structure [35] have been shown in a mathematically rigorous way to be exact approximations of the stochastic SIR dynamics in the limit of large node size. Our analytic PIDE models for SIR spread on networks with spatial structure were derived using similar assumptions, and the close match we obtain between stochastic and PIDE simulations as, e.g., in Fig. 5, indicates that there is a mathematical conjecture that our PIDE models are also exact in the limit of large node density. It may be possible to prove this mathematically using techniques similar to the ones used in [37, 38]. Furthermore, the RSN model defined by (1) is versatile as it allows flexible degree distributions and distance kernels. There are some further generalizations that were not considered in this article but offer compelling prospects for building realistic networks models that remain analysable. An obvious generalization is to allow $$\rho(\mathbf{x})$$ to describe an inhomogeneous spatial density of nodes. This could, for example, be used to model different population densities in cities and rural areas in the context of realistic spatial disease spreading models, or different densities of neurons in different parts of the brain. Similarly, instead of letting the distance kernel depend on nodal distances $$d_{uv}$$, one can consider more general kernels $$f(\mathbf{x}_u,\mathbf{x}_v)$$ that directly depend on the coordinates of the nodes. For example, this could model people living in cities preferentially connecting to people in the same and other cities, while connections within and between rural individuals could follow different patterns. This type of modelling is especially compelling in the era of big data where real data may be used to build analysable spatial random network models that faithfully mirror real-world spatial networks. For example, [10] studies the interplay between short-scale commuting flows and long-range airline traffic in shaping the spatiotemporal pattern of a global epidemic, and [11] found that during the 2013–2016 Ebola outbreak viral lineages moved according to a classic ‘gravity’ model, with more intense migration between larger and more proximate population centres. Using empirical data or inferred spatial connectivities, all these kinds of effects can be incorporated in our RSN models of (1), with clear potential for highly realistic models of disease propagation on spatial networks that retain the analytic tractability of the approach. Funding Joel C. Miller was funded by the Global Good Fund through the Institute for Disease Modeling and by a Larkins Fellowship from Monash University. Jamieson L. Kaiser was funded by an Australian Mathematical Sciences Institute Vacation Research Scholarship. John C. Lang acknowledges funding from the Natural Sciences and Engineering Research Council of Canada. Appendix A. Fast network generation At first sight, generating networks from the RSN class defined by Eq. (1) appears to be an $${\mathcal{O}}(N^2)$$ operation as each of the $$\binom{N}{2}$$ edges exists independently of the others. However, by modifying methods developed to generate Erdős–Rényi and Chung–Lu networks [44, 78] in linear time, we arrive at a much more efficient process. This makes large RSNs practical for simulation and analysis. We briefly outline the method, under the assumption that the distance kernel $$f$$ is bounded above and is monotonically decreasing. We divide the space $$V$$ into regions and order the nodes in each region by decreasing expected degree. We consider a region $$X_i$$, and define $$u$$ to be the first node in that region. The probability that $$u$$ will have an edge with any subsequent node is bounded above by $$q_0=\min(1,\kappa_u^2 f_{\text{max}}/{\left\langle {\kappa} \right \rangle})$$ where $$f_{\text{max}}= f(0)$$ is the maximum of $$f$$. We then choose a number $$r$$ from a geometric distribution with success probability $$q_0$$. So $$r$$ is chosen from a negative binomial distribution giving the number of unsuccessful trials (each having success probability $$q_0$$) before the first success. We set $$v_1$$ to be the $$r$$th node following $$u$$ in $$X_i$$. This is equivalent to $$v_1$$ being the first node chosen when sequentially considering each node with probability $$q_0$$. Thus the ‘candidate neighbour’ $$v_1$$ was chosen with probability $$q_0$$ which is at least $$p_{uv_1}=\min(1,\kappa_u \kappa_{v_1} f(d_{uv})/{\left\langle {\kappa} \right \rangle})$$. It is possible that $$v_1$$ is a “false positive”. To decide this, we generate a new random number uniformly from $$(0,1)$$, and if it is less than $$p_{uv_1}/q_0$$, we decide that $$v_1$$ was correctly chosen, and we add an edge between them. Otherwise $$v_1$$ is a false positive and we do not add the edge. We then enter an iterative step. After processing $$v_i$$, we set $$q_i = \min (1,\kappa_u \kappa_{v_i} f_{\text{max}}/{\left\langle {\kappa} \right \rangle})$$ and jump to the next candidate neighbour $$v_{i+1}$$. Again $$v_{i+1}$$ may be a false positive, and we place an edge between $$u$$ and $$v_{i+1}$$ with probability $$q_{uv_{i+1}}/q_i$$ where $$q_{uv_{i+1}} = \min (1,\kappa_u \kappa_{v_{i+1}} f(d_{uv_{i+1}})/{\left\langle {\kappa} \right \rangle})$$. As the iteration progresses, $$\kappa_{v_i}$$ decreases so $$q_i$$ decreases, and the jumps between successive $$v$$s become larger. The edges within each region are assigned in linear time. We next place edges between the node $$u$$ and other regions $$X_j$$. We find the minimum distance between $$u$$ and $$X_j$$ and use it to define $$f_{\text{max}}$$. We take $$w$$ to be the first node in $$X_j$$. We define $$q_0 = \min \{ 1, \kappa_u \kappa_w f_{\text{max}}/{\left\langle {\kappa} \right \rangle}\}$$. We choose $$v_1$$ from $$X_j$$ as before, and iterate. For $$X_j$$ farther from $$u$$ the jumps are larger. Appendix B. Relation to the Fisher–KPP equation We note that, under approximations that are appropriate in the case of a localized spatial kernel with $$\Theta$$ varying slowly in space we can convert (2) into the Fisher–KPP equation $$u_t = ru(1-u/K) + D u_{xx}$$. We demonstrate this in the 1D case. We start by assuming that $$S(x,0)$$ is approximately $$1$$. We then assume that $$f$$ is sufficiently localized and $$\Theta$$ varies sufficiently slowly that we can expand $$\Psi'(\Theta(x,t))$$ as a Taylor Series in $$x$$ to third order.   \begin{align*} \int_{-\infty}^\infty \Psi'(\Theta(\hat{x},t)) f(|\hat{x}-x|) \, \mathrm{d}\hat{x} = & \int_{-\infty}^\infty \Psi'(\Theta(x,t)) \, \mathrm{d}\hat{x}\\ & + \int_{-\infty}^\infty (\hat{x}-x)\frac{\partial {}}{\partial {x}}\Psi'(\Theta(x,t)) f(|\hat{x}-x|)\, \mathrm{d}\hat{x} \\ &+ \int_{-\infty}^\infty \frac{(\hat{x}-x)^2}{2}\left(\frac{\partial^{2} {}}{\partial x^{2}} \Psi'(\Theta(x,t)\right) f(|\hat{x}-x|) \, \mathrm{d}\hat{x} \\ &+ \cdots\\ \approx & \Psi'(\Theta(x,t)) + C\frac{\partial^{2} {}}{\partial x^{2}}\Psi'(\Theta(x,t)) \end{align*} where $$C = \int_{-\infty}^\infty y^2f(|y|)/2 \, \mathrm{d}y$$, and we have used the fact that $$\int_{-\infty}^\infty y f(|y|) \, \mathrm{d}y=0$$ by symmetry. So taking $$S(x,0)=1$$, our equation for $$\dot{\Theta}$$ is   \[ \dot{\Theta} = - \beta \Theta + \gamma(1-\Theta) + \beta \frac{\Psi'(\Theta)}{{\left\langle {\kappa} \right \rangle}} + \frac{\beta C}{{\left\langle {\kappa} \right \rangle}}\frac{\partial^{2} {}}{\partial x^{2}} \Psi'(\Theta) \] Now setting $$u = 1-\Theta$$ and assuming $$u$$ is small, we use further Taylor expansions and the fact that $$\Psi^{(n)}(1) = {\left\langle {\kappa^n} \right \rangle}$$ to find   \begin{align*} \dot{u} &= \beta (1-u) - \gamma u - \beta \frac{\Psi'(1-u)}{{\left\langle {\kappa} \right \rangle}} - \frac{\beta C}{{\left\langle {\kappa} \right \rangle}}\frac{\partial^{2} {}}{\partial x^{2}} \Psi'(1-u)\\ &\approx \beta (1-u) - \gamma u - \beta + \beta\frac{u{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}} - \beta \frac{u^2 {\left\langle {\kappa^3} \right \rangle}}{2{\left\langle {\kappa} \right \rangle}} - \frac{\beta C}{{\left\langle {\kappa} \right \rangle}}\frac{\partial^{2} {}}{\partial x^{2}}\Psi'(1-u)\\ &\approx \left(\beta \frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}-\beta-\gamma \right) u - \beta \frac{{\left\langle {\kappa^3} \right \rangle}}{2{\left\langle {\kappa} \right \rangle}}u^2 - \frac{\beta C}{{\left\langle {\kappa} \right \rangle}}\frac{\partial^{2} {}}{\partial x^{2}} \left({\left\langle {\kappa} \right \rangle} - u {\left\langle {\kappa^2} \right \rangle}\right)\\ &= r u(1-u/K) + D\frac{\partial^{2} {}}{\partial x^{2}} u \end{align*} for appropriately chosen $$r$$, $$K$$, and $$D$$. So the Fisher–KPP equation arises out of an approximation of our PIDE on RSNs. In general these approximations are not accurate when, e.g., $$\Theta$$ is not close to $$1$$, the variation in $$\Theta$$ is fast enough that the Taylor Series expansions are poor approximations, or the spatial kernel is not localized. However, this relation does suggest that behaviours found for the Fisher–KPP equation are likely to occur for disease spread in our RSNs as well. Footnotes 1 Technically the number of neighbours is binomially distributed with parameters $$N-1$$ and $$p=\kappa_u/(N-1)$$ which can be treated as Poisson-distributed with mean $$\kappa_u$$ for large $$N$$, or if we take $$\rho$$ as fixed and place nodes as a true Poisson process, then $$N$$ is a random variable, and the number of neighbours truly is Poisson distributed. References 1. Grenfell B. T., Bjørnstad O. N. & Kappey J. ( 2001) Travelling waves and spatial hierarchies in measles epidemics. Nature , 414, 716– 723. Google Scholar CrossRef Search ADS PubMed  2. Childs J. E., Curns A. T., Dey M. E., Real L. A., Feinstein l., Bjørnstad O. N. & Krebs J. W. ( 2000) Predicting the local dynamics of epizootic rabies among raccoons in the United States. Proc. Natl. Acad. Sci. U.S.A , 97, 13666– 13671. Google Scholar CrossRef Search ADS PubMed  3. Biek R., Henderson J. C., Waller L. A., Rupprecht C. E. & Real L. A. ( 2007) A high-resolution genetic signature of demographic and spatial expansion in epizootic rabies virus. Proc. Natl. Acad. Sci. U.S.A , 104, 7993– 7998. Google Scholar CrossRef Search ADS PubMed  4. Steck F. & Wandeler A. ( 1980) The epidemiology of fox rabies in europe. Epidemiol. Rev. , 2, 71– 96. Google Scholar CrossRef Search ADS PubMed  5. Baer G. M. ( 1991) The Natural History of Rabies. 2nd edn, Chapter 13-titled as ’fox rabies’  ( Baer G. M. ed.). Boca Raton, Ann Arbor and Boston: CRC Press. 6. Benavides J. A., Valderrama W. & Streicker D. G. ( 2016) Spatial expansions and travelling waves of rabies in vampire bats. In Proc. R. Soc. B , vol. 283, UK: The Royal Society, pp. 20160328. Google Scholar CrossRef Search ADS   7. Bos K. I., Schuenemann V. J., Golding G. B., Burbano H. A., Waglechner N., Coombes B. K., McPhee J. B., DeWitte S. N., Meyer M., Schmedes S. et al.   ( 2011) A draft genome of yersinia pestis from victims of the black death. Nature , 478, 506– 510. Google Scholar CrossRef Search ADS PubMed  8. Christakos G., Olea R. A. & Yu H.-L. ( 2007) Recent results on the spatiotemporal modelling and comparative analysis of black death and bubonic plague epidemics. Public Health , 121, 700– 720. Google Scholar CrossRef Search ADS PubMed  9. Marvel S. A., Martin T., Doering C. R., Lusseau D. & Newman M. E. J. ( 2013) The small-world effect is a modern phenomenon. Preprint arXiv:1310.2636. 10. Balcan D., Colizza V., Gonalves B., Hu H., Ramasco J. J. & Vespignani A. ( 2009) Multiscale mobility networks and the spatial spreading of infectious diseases. Proc. Natl. Acad. Sci. U.S.A , 106, 21484– 21489. Google Scholar CrossRef Search ADS PubMed  11. Dudas G., Carvalho L. M., Bedford T., Tatem A. J., Baele G., Faria N. R., Park D. J., Ladner J. T., Arias A., Asogun D., et al.   ( 2017) Virus genomes reveal factors that spread and sustained the ebola epidemic. Nature , 544, 309– 315. Google Scholar CrossRef Search ADS PubMed  12. Fine P. E. M. & Carneiro I. A. M. ( 1999) Transmissibility and persistence of oral polio vaccine viruses: implications for the global poliomyelitis eradication initiative. Am. J. Epidemiol. , 150, 1001– 1021. Google Scholar CrossRef Search ADS PubMed  13. Kroiss S. J., Famulare M., Lyons H., McCarthy K. A., Mercer L. D. & Chabot-Couture G. ( 2017) Evaluating cessation of the type 2 oral polio vaccine by modeling pre-and post-cessation detection rates. Vaccine , 35, 5674– 5681. Google Scholar CrossRef Search ADS PubMed  14. Newman M. E. J. ( 2002) Spread of epidemic disease on networks. Phys. Rev. E , 66, 016128. Google Scholar CrossRef Search ADS   15. Meyers L. A. ( 2007) Contact network epidemiology: Bond percolation applied to infectious disease prediction and control. Bull. Am. Math. Soc. , 44, 63– 86. Google Scholar CrossRef Search ADS   16. Pastor-Satorras R., Castellano C., Van Mieghem P. & Vespignani A. ( 2015) Epidemic processes in complex networks. Rev. Mod. Phys. , 87, 925. Google Scholar CrossRef Search ADS   17. Bullmore E. & Sporns O. ( 2009) Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. , 10, 186– 198. Google Scholar CrossRef Search ADS PubMed  18. O’Dea R., Crofts J. J. & Kaiser M. ( 2013) Spreading dynamics on spatially constrained complex brain networks. J. R. Soc. Interface , 10, 20130016. Google Scholar CrossRef Search ADS PubMed  19. Ercsey-Ravasz M., Markov N. T., Lamy C., Van Essen D. C., Knoblauch K., Toroczkai Z. & Kennedy H. ( 2013) A predictive network model of cerebral cortical connectivity based on a distance rule. Neuron , 80, 184– 197. Google Scholar CrossRef Search ADS PubMed  20. Horvát S., Gămănu R., Ercsey-Ravasz M., Magrou L., Gămănu B., Van Essen D. C., Burkhalter A., Knoblauch K., Toroczkai Z. & Kennedy H. ( 2016) Spatial embedding and wiring cost constrain the functional layout of the cortical network of rodents and primates. PLoS Biol. , 14, e1002512. Google Scholar CrossRef Search ADS PubMed  21. Barthélemy M. ( 2011) Spatial networks. Phys. Rep. , 499, 1– 101. Google Scholar CrossRef Search ADS   22. Auffinger A., Damron M. & Hanson J. ( 2017) 50 years of first passage percolation. Providence, RI: American Mathematical Society. 23. Chung F. & Lu L. ( 2002) Connected components in random graphs with given expected degree sequences. Ann. Comb. , 6, 125– 145. Google Scholar CrossRef Search ADS   24. Newman M. E. J. ( 2003) The structure and function of complex networks. SIAM Rev. , 45, 167– 256. Google Scholar CrossRef Search ADS   25. Centola D., Eguíluz V. M. & Macy M. W. ( 2007) Cascade dynamics of complex propagation. Phys. A Stat. Mech. Appl. , 374, 449– 456. Google Scholar CrossRef Search ADS   26. Janssen H.-K., Müller M. & Stenull O. ( 2004) Generalized epidemic process and tricritical dynamic percolation. Phys. Rev. E , 70, 026114. Google Scholar CrossRef Search ADS   27. Watts D. J. ( 2002) A simple model of global cascades on random networks. Proc. Natl. Acad. Sci. U.S.A , 99, 5766– 5771. Google Scholar CrossRef Search ADS PubMed  28. Miller J. C. ( 2016) Complex contagions and hybrid phase transitions. J. Complex Netw. , cnv021. 29. Melnik S., Hackett A., Porter M. A., Mucha P. J. & Gleeson J. P. ( 2011) The unreasonable effectiveness of tree-based theory for networks with clustering. physical review E , 83, 036112. Google Scholar CrossRef Search ADS   30. Kermack W. O. & McKendrick A. G. ( 1927) A contribution to the mathematical theory of epidemics. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences , vol. 115. The Royal Society, pp. 700– 721. 31. Anderson R. M. & May R. M. ( 1991) Infectious Diseases of Humans . Oxford: Oxford University Press. 32. Kiss I. Z., Miller J. C. & Simon P. L. ( 2017) Mathematics of Epidemics on Networks: From Exact to Approximate Models . IAM. New York: Springer. Google Scholar CrossRef Search ADS   33. Abuelezam N. N., Rough K. & Seage III G. R. ( 2013) Individual-based simulation models of HIV transmission: Reporting quality and recommendations. PLoS One , 8, e75624. Google Scholar CrossRef Search ADS PubMed  34. Eaton J. W., Johnson L. F., Salomon J. A., Bärnighausen T., Eran B., Bershteyn A., Bloom D. E., Cambiano V., Fraser C., Hontelez J. A. C., et al.   ( 2012) HIV treatment as prevention: systematic comparison of mathematical models of the potential impact of antiretroviral therapy on HIV incidence South Africa. PLoS Med. , 9, e1001245. Google Scholar CrossRef Search ADS PubMed  35. Miller J. C., Slim A. C. & Volz E. M. ( 2012) Edge-based compartmental modelling for infectious disease spread. J. R. Soc. Interface , 9, 890– 906. Google Scholar CrossRef Search ADS PubMed  36. Molloy M. & Reed B. ( 1995) A critical point for random graphs with a given degree sequence. Random Struct. & Algorithms , 6, 161– 180. Google Scholar CrossRef Search ADS   37. Janson S., Luczak M. & Windridge P. ( 2014) Law of large numbers for the SIR epidemic on a random graph with given degrees. Random Struct. & Algorithms , 45, 724– 761. Google Scholar CrossRef Search ADS   38. Decreusefond L., Dhersin J.-S., Moyal P. & Tran V. C. ( 2012) Large graph limit for an SIR process in random network with heterogeneous connectivity. Ann. Appl. Prob. , 22, 541– 575. Google Scholar CrossRef Search ADS   39. Diekmann O. ( 1979) Run for your life. A note on the asymptotic speed of propagation of an epidemic. J. Differ. Equ. , 33, 58– 73. Google Scholar CrossRef Search ADS   40. Diekmann O. ( 1978) Thresholds and travelling waves for the geographical spread of infection. J. Math. Biol. , 6, 109– 130. Google Scholar CrossRef Search ADS PubMed  41. Mollison D. ( 1972) The rate of spatial propagation of simple epidemics. In Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 3: Probability Theory . Berkeley and Los Angeles: The Regents of the University of California. 42. Bollobás B., Janson S. & Riordan O. ( 2007) The phase transition in inhomogeneous random graphs. Random Struct. & Algorithms , 31, 3– 122. Google Scholar CrossRef Search ADS   43. Norros I. & Reittu H. ( 2006) On a conditionally Poissonian graph process. Adv. Appl. Prob. , 38, 59– 75. Google Scholar CrossRef Search ADS   44. Miller J. C. & Hagberg A. ( 2011) Efficient generation of networks with given expected degrees. Proceedings of the 8th International Workshop on Algorithms and Models for the Web Graph . Berlin: Springer, pp. 115– 126. Google Scholar CrossRef Search ADS   45. Söderberg B. ( 2002) General formalism for inhomogeneous random graphs. Physical Review E , 66, 066121. Google Scholar CrossRef Search ADS   46. Gilbert E. N. ( 1959) Random graphs. Ann. Math. Stat. , 30, 1141– 1144. Google Scholar CrossRef Search ADS   47. Watts D. J. & Strogatz S. H. ( 1998) Collective dynamics of small-world networks. Nature , 393, 440– 442. Google Scholar CrossRef Search ADS PubMed  48. Newman M. E. J. & Watts D. J. ( 1999) Renormalization group analysis of the small-world network model. Phy. Lett. A , 263, 341– 346. Google Scholar CrossRef Search ADS   49. Kleinberg J. M. ( 2002) Small-world phenomena and the dynamics of information. In Advances in neural information processing systems ., vol. 15. Boston: MIT Press, pp. 431– 438. 50. Waxman B. M. ( 1988) Routing of multipoint connections. IEEE J. Selected Areas Commun. , 6, 1617– 1622. Google Scholar CrossRef Search ADS   51. Newman C. M. & Schulman L. S. ( 1986) One dimensional $$1/| j- i|^s$$ percolation models: The existence of a transition for $$s \leq 2$$. Commun. Math. Phys. , 104, 547– 571. Google Scholar CrossRef Search ADS   52. Biskup M. ( 2004) On the scaling of the chemical distance in long-range percolation models. Ann. Prob. , 32, 2938– 2977. Google Scholar CrossRef Search ADS   53. Trapman P. ( 2010) The growth of the infinite long-range percolation cluster. Ann. Prob. , 38, 1583– 1608. Google Scholar CrossRef Search ADS   54. Deijfen M., van der Hofstad R. & Hooghiemstra G. ( 2013) Scale-free percolation. In Annales de l’Institut Henri Poincaré, Probabilités et Statistiques , vol. 49. Institut Henri Poincaré, pp. 817– 838. Google Scholar CrossRef Search ADS   55. Deprez P. & Wüthrich M. V. ( 2015) Networks, random graphs and percolation. In Theoretical Aspects of Spatial-Temporal Modeling . Tokyo: Springer, pp. 95– 124. Google Scholar CrossRef Search ADS   56. Vladimirov N., Traub R. D. & Tu Y. ( 2011) Wave speed in excitable random networks with spatially constrained connections. PLoS One , 6, e20536. Google Scholar CrossRef Search ADS PubMed  57. Barnett L., Di Paolo E. & Bullock S. ( 2007) Spatially embedded random networks. Phys. Rev. E , 76, 056115. Google Scholar CrossRef Search ADS   58. Kosmidis K., Havlin S. & Bunde A. ( 2008) Structural properties of spatially embedded networks. Europhys. Lett. , 82, 48005. Google Scholar CrossRef Search ADS   59. Robins G., Pattison P., Kalish Y. & Lusher D. ( 2007) An introduction to exponential random graph ($$p*$$) models for social networks. Soc. Netw. , 29, 173– 191. Google Scholar CrossRef Search ADS   60. Wong L. H., Pattison P. & Robins G. ( 2006) A spatial model for social networks. Phys. A Stat. Mech. Appl. , 360, 99– 120. Google Scholar CrossRef Search ADS   61. Haenggi M., Andrews J. G., Baccelli F., Dousse O. & Franceschetti M. ( 2009) Stochastic geometry and random graphs for the analysis and design of wireless networks. IEEE J. Selected Areas Commun. , 27, 1029– 1046. Google Scholar CrossRef Search ADS   62. Brockmann D. & Helbing D. ( 2013) The hidden geometry of complex, network-driven contagion phenomena. Science , 342, 1337– 1342. Google Scholar CrossRef Search ADS PubMed  63. Hoff P. D., Raftery A. E. & Handcock M. S. ( 2002) Latent space approaches to social network analysis. J. Am. Stat. Assoc. , 97, 1090– 1098. Google Scholar CrossRef Search ADS   64. Matias C. & Robin S. ( 2014) Modeling heterogeneity in random graphs through latent space models: a selective review. ESAIM Proc. Surv. , 47, 55– 74. Google Scholar CrossRef Search ADS   65. Serrano M. A., Krioukov D. & Boguná M. ( 2008) Self-similarity of complex networks and hidden metric spaces. Phys. Rev. Lett. , 100, 078701. Google Scholar CrossRef Search ADS PubMed  66. Volz E. M. ( 2008) SIR dynamics in random networks with heterogeneous connectivity. J. Math. Biol. , 56, 293– 310. Google Scholar CrossRef Search ADS PubMed  67. Miller J. C. ( 2014) Epidemics on networks with large initial conditions or changing structure. PLoS One , 9, e101421. Google Scholar CrossRef Search ADS PubMed  68. Fisher R. A. ( 1937) The wave of advance of advantageous genes. Ann. Eugenics , 7, 355– 369. Google Scholar CrossRef Search ADS   69. Kolmogorov A. N., Petrovskii I. G. & Piskunov N. S. ( 1991) A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem. In: Tikhomirov V. M. ed. Selected Works of A. N. Kolmogorov . Dordrecht: Springer. 70. Coville J. & Dupaigne L. ( 2005) Propagation speed of travelling fronts in non local reaction–diffusion equations. Nonlinear Anal. Theory Meth. Appl. , 60, 797– 819. Google Scholar CrossRef Search ADS   71. Berestycki H., Nadin G., Perthame B. & Ryzhik L. ( 2009) The non-local Fisher-KPP equation: travelling waves and steady states. Nonlinearity , 22, 2813. Google Scholar CrossRef Search ADS   72. Hallatschek O. & Fisher D. S. ( 2014) Acceleration of evolutionary spread by long-range dispersal. Proc. Natl. Acad. Sci. U.S.A , 111, E4911– E4919. Google Scholar CrossRef Search ADS PubMed  73. Miller J. C. & Hagberg A. ( 2011) Efficient generation of networks with given expected degrees. In Proceedings of the 8th International Workshop on Algorithms and Models for the Web Graph . pp. 115– 126. 74. Panja D. ( 2004) Effects of fluctuations on propagating fronts. Phys. Rep. , 393, 87– 174. Google Scholar CrossRef Search ADS   75. Penrose M. ( 2003) Random Geometric Graphs , vol. 5. Oxford: Oxford University Press. Google Scholar CrossRef Search ADS   76. Joos F., Perarnau G., Rautenbach D. & Reed B. ( 2016) How to determine if a random graph with a fixed degree sequence has a giant component. IEEE 57th Annual Symposium on Foundations of Computer Science , pp. 695– 703. 77. Bollobás B., Janson S. & Riordan O. ( 2007) Spread-out percolation in d. Random Struct. & Algorithms , 31, 239– 246. Google Scholar CrossRef Search ADS   78. Batagelj V. & Brandes U. ( 2005) Efficient generation of large random networks. Phys. Rev. E , 71, 036113. Google Scholar CrossRef Search ADS   © The authors 2018. Published by Oxford University Press. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Journal of Complex Networks Oxford University Press

Analytic models for SIR disease spread on random spatial networks

Loading next page...
 
/lp/ou_press/analytic-models-for-sir-disease-spread-on-random-spatial-networks-WYER5LFcLi
Publisher
Oxford University Press
Copyright
© The authors 2018. Published by Oxford University Press.
ISSN
2051-1310
eISSN
2051-1329
D.O.I.
10.1093/comnet/cny004
Publisher site
See Article on Publisher Site

Abstract

Abstract We study the propagation of susceptible-infectious-recovered (SIR) disease on random networks with spatial structure. We consider random spatial networks (RSNs) with heterogeneous degrees, where any two nodes are connected by an edge with a probability that depends on their spatial distance and expected degrees. This is a natural spatial extension of the well-known Erdős–Rényi and Chung–Lu random graphs. In the limit of high node density, we derive partial integro-differential equations governing the spread of SIR disease through an RSN. We demonstrate that these nonlocal evolution equations quantitatively predict stochastic simulations. We analytically show the existence of travelling wave solutions on the RSNs. If the distance kernel governing edge placement in the RSNs decays slower than exponential, the population-scale dynamics are dominated by long-range hops followed by local spread of travelling waves. This provides a theoretical modelling framework for disease propagation on networks with spatial structure, extending analytic SIR modelling capabilities from networks without spatial structure to the important real-life case of networks with spatial structure. Potential applications of this modelling framework include studying the interplay between long-range and short-range transmissions in designing polio vaccination strategies, wavelike spread of the plague in medieval Europe, and recent observations of how epidemics like Ebola evolve in modern connected societies. 1. Introduction The spread of infectious disease through human and animal populations exhibits a range of patterns in space and time. In the pre-vaccination era, Measles in the UK spread out from London in a coherent spatial pattern [1]. Other diseases exhibiting coherent spatial expansion include rabies in raccoons [2, 3], foxes [4, 5] and vampire bats [6] as well as the Black Death of 1347–1351 in Europe, which claimed an estimated 30–50% of the European population [7–9]. These coherent spatial patterns seem the norm in many epizootics as well as historical human epidemics. However, more recent human epidemics exhibit different patterns. The long-range connections in modern societies play increasingly important roles in disease propagation [10]. These lead to the appearance of multiple spreading hotspots of disease. The 2013–2016 epidemic of Ebola in West Africa had significant long-range hops: sequencing of over 1600 Ebola virus genomes reveals a heterogeneous and spatially dissociated collection of transmission clusters of varying size, duration and connectivity [11] (Video 1 of [11] illustrates such a pattern of spread in the Ebola epidemic that we may refer to as ‘hop-and-spread’ dynamics). Similarly, SARS and pandemic 2009 H1N1 influenza showed significant local spread, but the global dynamics were dominated by sporadic long-range transmission events. A further example is the influence of spatial structure on polio eradication. Oral polio vaccine (OPV) is a live, weakly transmitted virus which—such as the wild virus—can cause paralysis and can evolve to be more transmissible [12]. Wild type-2 polio is extinct, so ideally we would like to remove OPV-2 vaccination, but cessation comes with the risk that a transmission chain of the vaccine strain may persist until enough unvaccinated children are born to allow it to cause long-scale epidemics [13]. Routine OPV-2 vaccination has stopped, and the vaccine is currently being used in response to identified circulation of vaccine-derived strains. This may eliminate the circulating strain, but with a risk of seeding new transmission chains. Understanding the interplay between long-range and short-range transmissions is important in designing strategies for localized OPV-2 use, and in informing decisions about stockpile maintenance. Network models are a valuable tool for investigating how a population’s contact structure governs the spread of infectious disease [14–16]. They also appear in many other contexts, such as understanding activation patterns in neurons. In many networks of interest, there is an underlying geometry to the location of nodes, and nearby nodes are more likely to be connected than nodes that are separated by a large distance [17–21]. Outside of first-passage percolation [22] where the networks are generally lattices, there is relatively little theoretical study of how a network’s spatial structure affects spreading processes. Commonly used spatial network models have a number of weaknesses. In particular, they are not amenable to analytic investigation, and many of them have subtle preferences for specific directions. By way of contrast, non-spatial network models such as Chung and Lu [23] and Configuration model networks [24] provide significant insight. This is largely because their structure has allowed for analytic exploration. They are ‘locally tree-like’: for fixed cycle length $$D$$, the probability that a random node is in a cycle of length $$D$$ or less goes to zero as the number of nodes grows. The dynamics of many spreading processes, for example, complex contagions [25], the generalized epidemic process [26] and the Watts threshold model [27] can be studied exactly in locally tree-like networks [28], and even when the tree-like hypothesis breaks down, spreading processes on networks are often accurately predicted by the unclustered assumptions [29]. To study these in networks with spatial structure, we would like a spatial network model that contains the locally tree-like behaviour as a limiting case. The susceptible-infectious-recovered (SIR) model is typically formulated as a system of ordinary differential equations (ODEs) for well-mixed populations [30, 31]. To understand how network effects such as social heterogeneity and nonzero partnership duration affect spreading dynamics, the SIR model is often simulated stochastically [32–34]. Recently, [35] introduced the edge-based compartmental modelling (EBCM) technique that leads to low-dimensional ODE systems that match the stochastic dynamics in the limit of large network size for random configuration model networks [24] (also known as Molloy–Reed graphs [36]). It has been shown mathematically in a rigorous way that the ODE approximation of [35] is indeed exact in the limit [37, 38]. The ODE models allow analytic investigation of the SIR disease processes on networks. The networks of [35], however, do not have spatial structure, and there are many important aspects of SIR disease propagation related to the geographical spread of infection that depend on the spatial aspects of population structure [39–41], in addition to the network structure. However, to date, studies of the spread of disease in spatial networks have relied on simulation because the network models used do not lead to a simple analytic model. In this article, we fill this gap, and extend the analytic edge-based compartmental modelling capabilities of [28, 35] for non-spatial networks to spreading processes on networks with spatial structure. We use a spatial network model that forms a natural extension of the well-known Erdős–Rényi and Chung–Lu random graphs, adding the crucial element of spatial structure while keeping the model amenable to mathematical analysis. This allows us to analytically study how spatial structure affects spreading dynamics. In particular, it enables studying the interplay between long-range and short-range transmissions of SIR disease on networks with spatial structure, and how this affects travelling waves and hop-and-spread dynamics, with applications to understanding propagation of diseases like the plague in medieval Europe and Ebola in modern connected societies, and designing polio vaccination strategies. The remainder of this article is structured as follows. In Section 2, we investigate analytic models for SIR disease propagation on random spatial networks (RSNs). We discuss the RSN class of random networks with spatial structure that we consider in this article, and derive an analytic partial integro-differential equation (PIDE) model to quantitatively describe SIR dynamics on RSNs in the limit of large node density. We show analytically that SIR disease propagation on RSNs in this limit allows for exact travelling wave solutions (i.e., nonlinear waves that retain their shape while propagating) characterised by the nonlocal evolution equations of PIDE type. We study the analytic properties of these travelling wave solutions and conditions under which they exist. In Section 3, we apply stochastic SIR simulation to examples of RSNs, and demonstrate that the nonlocal evolution equations quantitatively predict stochastic simulations. We identify travelling waves in the simulations and verify the analytic properties derived from the PIDE. The interplay between short-range and long-range connections reveals dynamics in which travelling waves emanate from focal points that result from long-range hops. Section 4 discusses the application of these new models to problems in disease propagation and concludes. 2. SIR disease propagation on RSNs: analytic models and travelling waves We first introduce the class of RSNs that we consider in this article. We then derive differential equations that govern the spatial dynamics of SIR disease on RSNs in the limit of large node density and demonstrate the existence of nonlinear travelling wave solutions that retain their shapes. We derive analytic properties of the resulting differential equations, including properties of the travelling waves. These models are used in Section 3 to study the dynamics of SIR disease on RSNs. 2.1 The RSN class To study disease propagation on spatial networks, we consider a class of RSNs modelled after Chung–Lu networks [23] (formally a special case of the inhomogeneous random graph class [42]). We place nodes uniformly at random into some region $$V$$ with density $$\rho$$. For our purposes $$\rho = N/|V|$$, where $$N$$ is the total number of nodes and $$|V|$$ is the (assumed finite) volume of the space $$V$$. Each node $$u$$ is independently assigned an expected degree $$\kappa_u$$ from some distribution $$P(\kappa)$$. The average of $$\kappa$$ is denoted $${\left\langle {\kappa} \right \rangle}$$. A $$u$$–$$v$$ edge is created with probability   \begin{equation} p_{uv} = \min\left(\kappa_u\kappa_v \frac{f(d_{uv})}{\rho {\left\langle {\kappa} \right \rangle}}, 1\right). \label{eq:p} \end{equation} (1) We assume that the distance kernel$$f$$ depends only on distance, is non-negative, and integrates to $$1$$: $$\int_V f(|\mathbf{x}-\mathbf{x}_0|) \, \mathrm{d}\mathbf{x} = 1$$ for all $$\mathbf{x}_0$$. Typically, we also assume that $$f$$ decreases monotonically with distance. An alternate function that would generally yield similar behaviour in the large $$\rho$$ limit would be $$p_{uv} = 1- e^{-\kappa_u \kappa_v f(d_{uv})/\rho {\left\langle {\kappa} \right \rangle}}$$, following [43]. If $$p_{uv}<1$$ always then the expected degree of a node $$u$$ is   \[ \int_V \int_0^\infty \rho f(|\mathbf{x}_v-\mathbf{x}_u|) \kappa_u \kappa_v \frac{P(\kappa_v)}{\rho {\left\langle {\kappa} \right \rangle}} \mathrm{d}\kappa_v \, \mathrm{d}\mathbf{x}_v = \kappa_u \, . \] In the large network limit, the degree of $$u$$ is Poisson distributed with mean $$\kappa_u$$. This model offers considerable flexibility: The distance kernel $$f$$ can be tuned to generate a wide range of spatial structure. The expected degree distribution allows us to tune the distribution of realized degrees. Many further generalizations (presented in the discussion in Section 4) emerge naturally. Figure 1 shows part of an RSN with $$10^4$$ nodes in $$[0,1]\times[0,1]$$ with an imposed bimodal degree distribution and a Gaussian distance kernel. Fig. 1. View largeDownload slide An example RSN and its properties. The distance kernel is a Gaussian, $$f(d) = \exp(-d^2/2\sigma^2)/2\pi\sigma^2$$ with $$\sigma = 0.03$$. The imposed distribution of expected degrees is $$P(2)=P(15)=0.5$$. The density is $$\rho=10000$$. One node and its neighbours are highlighted. A random network without spatial structure would exhibit neighbours throughout the domain. Fig. 1. View largeDownload slide An example RSN and its properties. The distance kernel is a Gaussian, $$f(d) = \exp(-d^2/2\sigma^2)/2\pi\sigma^2$$ with $$\sigma = 0.03$$. The imposed distribution of expected degrees is $$P(2)=P(15)=0.5$$. The density is $$\rho=10000$$. One node and its neighbours are highlighted. A random network without spatial structure would exhibit neighbours throughout the domain. A practical network model needs to be computationally feasible to generate. Appendix A describes an efficient generation algorithm for RSN networks, based on an $${\mathcal{O}}(N)$$ algorithm for Chung–Lu networks [44]. If the density is high enough that $$p_{uv}<1$$ for all node pairs, some important properties emerge. Consider a node $$u$$ at position $$\mathbf{x}$$ with a given $$\kappa_u$$, and consider network realizations containing $$u$$ but with different values of $$\rho$$. As $$\rho$$ increases, the distribution of the number of neighbours of $$u$$ in the network and the distribution of their locations does not change1. However, the probability that any two neighbours of $$u$$ are themselves neighbours of one another scales like $$1/\rho$$. Thus we anticipate that these networks will have similar spatial structure, but their clustering becomes negligible in the large $$\rho$$ limit. So the networks are locally tree-like at large $$\rho$$. This will allow our later analytic approaches to work. 2.1.1 Related network models We briefly review random network models that are related to RSNs. RSNs are a subclass of inhomogeneous random graphs [42, 45] and generalize many widely-studied models. They incorporate both degree heterogeneity and spatial structure in a way that, crucially, remains analytically tractable. They are a natural spatial extension of the well-known Erdős–Rényi and Chung–Lu random graphs. In a Chung–Lu network, an edge between $$u$$ and $$v$$ exists with probability $$\kappa_u\kappa_v/(N{\left\langle {\kappa} \right \rangle})$$ [23]. We recover this by setting $$f(d)=1/|V|$$ in Eq. (1). If we further set $$\kappa_u=\kappa$$ to be fixed for all $$u$$, we arrive at the “Erdős–Rényi” graph model, introduced by Gilbert [46]. RSNs are related to some other network models with spatial structure. There are a number of models which place nodes into a lattice and then assign edges by a rule that is at least in part determined by the distance between nodes. These include small-world networks [47, 48], navigable small-world networks [49], Waxman graphs [50], long-range percolation [51–53], scale-free percolation [54] (in particular [55] which uses a model very similar to our own) and spatially-constrained Erdős–Rényi graphs [56]. Most of the lattice-based models can be generated from the RSN model if we place the nodes on lattice points rather than uniformly at random. The widely-used random geometric graphs arise as a subclass of RSNs if we take uniform $$\kappa$$ and appropriately choose $$f$$. Furthermore, the class of spatially embedded random networks [57, 58], arises from a constant $$\kappa$$ and a power-law for $$f$$. The exponential random graph model does not have a clear connection to our model. It can handle degree heterogeneity [59] and some spatial structure [60], but typically the actual networks considered are small because they are expensive to generate. Many of these spatial networks are used to understand patterns emerging in the brain [18] or wireless sensor networks [61]. These models often have difficult-to-control correlations because of significant clustering. The microscopic structure makes analytic investigation of dynamic processes difficult. Due to this last weakness, almost all studies of spreading processes in networks with spatial structure are limited to simulation [62]. When equation-based analysis is attempted, it lacks quantitative agreement with simulations [56]. A major advantage of RSNs compared with other networks with spatial structure is their suitability to analytic results in the high-density limit. Although we demonstrate this for disease spread, many analytic techniques used to study other spreading processes in random networks will also apply to RSNs because of their locally tree-like structure in the high-density limit. We note however, that in some cases the high-density limit may not be achieved or significant clustering may exist, so the analytic model will not be accurate. In such cases, it may be possible to carry out an asymptotic expansion to identify corrections due to the lower density. We do not perform a detailed investigation of this. Throughout we assume that the nodes are embedded in a known geometry and a random network is generated for which the edges are assigned based on their distance and expected degrees. A closely related approach is given by latent space network models [63, 64]. In these models a network is given as the input, and the researcher attempts to find what underlying embedding geometry would most likely lead to such a network assuming a network generation algorithm similar to ours. There has been a significant amount of research based on this context, and many of these concepts also arise in studies on self-similar networks in hidden metric spaces [65]. Although much of the underlying mathematics in the network models is the same, our goal in using RSNs is to understand how processes would spread in a random network embedded in a given geometry, rather than to identify what geometry led to the creation of a given network. 2.2 Analytic governing equations for SIR disease on RSNs Our goal is to derive differential equations that accurately describe SIR disease dynamics on RSNs. We have noted that in the high-density limit RSNs have the crucial property of being locally tree-like, which allows us to apply and extend analytic tools that have been applied to other locally tree-like non-spatial networks like Chung–Lu or configuration model networks [35]. We derive analytic equations governing the deterministic dynamics (i.e., assuming that the number of nodes is sufficiently large for stochastic effects to become negligibly small). Our approach is of broad potential relevance, since we expect that analytic models on non-spatial graphs for processes other than SIR can be adapted to RSNs in a similar manner. In many random network classes, it is possible to develop powerful mathematical approaches by using their locally tree-like property: as $$N \to \infty$$ the clustering scales like $$1/N$$. This also occurs for our random spatial graphs: as $$\rho \to \infty$$, the clustering scales like $$1/\rho$$. We use this to develop analytic equations for SIR disease spread in the deterministic limit, following [35, 66, 67], but extended to the spatial setting. In an SIR epidemic, nodes begin susceptible and may be infected by infected neighbours (with rate $$\beta$$ per edge). Eventually infected nodes recover (with rate $$\gamma$$ per node). We assume that the population-scale dynamics are deterministic. That is, recognizing that the underlying process is stochastic, we assume that, as a function of $$\mathbf{x}$$, the proportion infected at some later time $$t$$ is uniquely determined from the initial proportion infected (as a function of $$\mathbf{x}$$). This assumption becomes reasonable in the limit of a large network, and is effectively the continuum assumption. Figure 2 illustrates that the spatio-temporal dynamics indeed become increasingly deterministic as the node density increases. Fig. 2. View largeDownload slide Disease spread in spatial networks in $$[0,1]\times[0.1]$$. The networks have $$P(3)=P(15)=1/2$$ and only short range connections [$$f(d)=1/0.01^2\pi$$ for $$d<0.01$$ and $$0$$ otherwise]. The value of $$\rho$$ varies, but has no effect on the bimodal degree distribution. Top: $$\rho=10^5$$, Middle: $$\rho=5\times 10^5$$, and Bottom: $$\rho=10^6$$. Disease is introduced to all nodes in a region of radius $$0.05$$ about the centre, with $$\beta=1/3$$ and $$\gamma=1$$. We plot the status of nodes at time $$t=12$$, with a detailed view of a subregion with the colours denoting a recovered node, an infected node or a susceptible node. As $$\rho$$ increases, the spread becomes more coherent and a deterministic limiting behaviour emerges. Fig. 2. View largeDownload slide Disease spread in spatial networks in $$[0,1]\times[0.1]$$. The networks have $$P(3)=P(15)=1/2$$ and only short range connections [$$f(d)=1/0.01^2\pi$$ for $$d<0.01$$ and $$0$$ otherwise]. The value of $$\rho$$ varies, but has no effect on the bimodal degree distribution. Top: $$\rho=10^5$$, Middle: $$\rho=5\times 10^5$$, and Bottom: $$\rho=10^6$$. Disease is introduced to all nodes in a region of radius $$0.05$$ about the centre, with $$\beta=1/3$$ and $$\gamma=1$$. We plot the status of nodes at time $$t=12$$, with a detailed view of a subregion with the colours denoting a recovered node, an infected node or a susceptible node. As $$\rho$$ increases, the spread becomes more coherent and a deterministic limiting behaviour emerges. We now derive analytic equations that faithfully describe the large $$\rho$$ SIR dynamics on RSNs. This will allow us to study and derive important properties of SIR disease propagation on networks with spatial structure. We start with the observation that the following four questions have identical answers in the large density limit, if indeed the population-scale dynamics are deterministic: As a function of location $$\mathbf{x}$$, what fraction of nodes in an infinitesimal region about $$\mathbf{x}$$ are susceptible, infected, or recovered at time $$t$$? What is the probability a random node is susceptible, infected, or recovered at time $$t$$, as a function of location $$\mathbf{x}$$? What is the probability a random node is susceptible, infected or recovered at time $$t$$, as a function of location $$\mathbf{x}$$, given the system’s state at time $$0$$? What is the probability a randomly chosen test node $$u$$ is susceptible, infected or recovered given the system’s state at time $$0$$ if we prevent $$u$$ from transmitting to its neighbours? The first two questions are clearly equivalent. The second and third are equivalent because we assume that the population-scale dynamics are deterministic. The third and fourth are equivalent because as long as $$u$$ is susceptible, the fact it does not transmit is irrelevant, and once $$u$$ is infected, its recovery time is not affected by any transmissions it causes. So the requirement that $$u$$ does not transmit does not influence the evolution of the state of $$u$$; it influences the states of neighbours of $$u$$, but only after it is too late for them to influence $$u$$. The process we analyse when preventing $$u$$ from transmitting (question 4) is different from the process in questions 1–3, but in question 4, we are not asking about the statuses of all nodes in the network. Rather we only ask about the status of node $$u$$, and for both processes, this is the same. Note that the equivalence of questions 3 and 4 holds specifically for SIR disease; it is crucial for our argument that once an individual begins to alter the status of its neighbours, its own future status can never be affected by those neighbours. We address question 4 using probabilistic tools. We define a test node $$u$$ which is randomly selected in the network at location $$\mathbf{x}$$ and prevented from transmitting to its neighbours. Under the assumption that $$\rho$$ is large, we can assume independence of $$u$$’s neighbours because clustering is negligible and because $$u$$ does not form a path of infection between its neighbours, precisely because we prevent it from transmitting. We seek to find the probability $$u$$ is susceptible, from which we deduce the probability it is infected or recovered. We pass to a continuum limit and define $$S(\mathbf{x},t)$$, $$I(\mathbf{x},t)$$ and $$R(\mathbf{x},t)$$ to be the probability that a test node $$u$$ placed at $$\mathbf{x}$$ would be susceptible, infected or recovered at time $$t$$. We similarly define $$\Phi_S(\mathbf{x},t)$$, $$\Phi_I(\mathbf{x},t)$$, $$\Phi_R(\mathbf{x},t)$$, and $$\Theta(\mathbf{x},t)$$ such that: $$\Phi_S$$ is the probability a random neighbour of the test node $$u$$ is susceptible at time $$t$$, $$\Phi_I$$ is the probability the neighbour is infected but has not transmitted to $$u$$, $$\Phi_R$$ is the probability the neighbour has recovered without transmitting to $$u$$ and $$\Theta=\Phi_S+\Phi_I+\Phi_R$$ is, thus, the probability that a random neighbour of node $$u$$ at position $$\mathbf{x}$$ has not transmitted infection to $$u$$ by time $$t$$. We assume infection is introduced to the population at time $$t=0$$ with $$S(\mathbf{x},0)$$ denoting the probability a node at $$\mathbf{x}$$ would be susceptible. We assume all other nodes are initially susceptible $$S(\mathbf{x},0) = 1- I(\mathbf{x},0)$$. Because $$k$$, the number of neighbours $$u$$ has, is Poisson distributed with mean $$\kappa_u$$ (the probability of a given $$k$$ is $$e^{-\kappa_u}\kappa_u^k/k!$$), the probability $$u$$ is susceptible at a later time $$t$$ is $$S(\mathbf{x},0) \sum_{k=0}^\infty e^{-\kappa_u}\kappa_u^k \Theta^k/k! = S(\mathbf{x},0) \exp[-\kappa_u(1-\Theta)]$$. Note that we have used here that the neighbours are independent (since we prevent $$u$$ from transmitting). We do not know the value of $$\kappa_u$$, so we must average over $$\kappa$$. The probability $$u$$ is susceptible is   \begin{align*} S &= S(\mathbf{x},0) \, \Psi(\Theta(\mathbf{x},t))\\ & = S(\mathbf{x},0) \int_0^\infty e^{-\kappa(1-\Theta(\mathbf{x},t))} P(\kappa) \, \mathrm{d}\kappa, \end{align*} with $$\Psi(\Theta(\mathbf{x},t)) \equiv \int_0^\infty e^{-\kappa(1-\Theta(\mathbf{x},t))} P(\kappa) \, \mathrm{d}\kappa$$. We augment this with $$I = 1-S-R$$ and $$\dot{R} = \gamma I$$. We must still find $$\Theta(\mathbf{x},t)$$. We turn to the flow diagram in Fig. 3. Once a neighbour $$v$$ of the test node $$u$$ becomes infected, whether or not $$v$$ transmits infection to $$u$$ and if so, at what time the transmission occurs, is independent of anything else that happens in the population. Thus we can immediately calculate that the flux from $$\Phi_I$$ to $$1-\Theta$$ is $$\beta\Phi_I$$, so $$\dot{\Theta} = -\beta \Phi_I$$. Similarly we find that the flux to $$\Phi_R$$ is $$\gamma \Phi_I$$. Since the flux to $$1-\Theta$$ and $$\Phi_R$$ are both proportional to $$\Phi_I$$ and both are $$0$$ at $$t=0$$, we can conclude that $$\Phi_R = \gamma(1-\Theta)/\beta$$. However, the rate at which $$v$$ becomes infected is more difficult. It will be easier to directly calculate $$\Phi_S$$ in terms of $$\Theta$$ and use $$\Phi_I = \Theta-\Phi_S-\Phi_R$$ to give an equation for $$\dot{\Theta}$$ in terms of $$\Theta$$. Fig. 3. View largeDownload slide A flow diagram leading to the governing equation for $$\Theta(\mathbf{x},t)$$. The compartments $$\Phi_S$$, $$\Phi_I$$ and $$\Phi_R$$ make up $$\Theta$$ and represent the probability that a random partner has not transmitted to $$u$$ and is susceptible, infected or recovered, respectively. Fig. 3. View largeDownload slide A flow diagram leading to the governing equation for $$\Theta(\mathbf{x},t)$$. The compartments $$\Phi_S$$, $$\Phi_I$$ and $$\Phi_R$$ make up $$\Theta$$ and represent the probability that a random partner has not transmitted to $$u$$ and is susceptible, infected or recovered, respectively. To find $$\Phi_S$$, we consider the possible neighbours of $$u$$ and calculate their probability of being susceptible. The probability a node $$v$$ at $$\hat{\mathbf{x}}$$ is a neighbour of $$u$$ is proportional to $$\kappa_v$$ and to $$f(|\hat{\mathbf{x}}-\mathbf{x}|)$$. In turn, the probability $$v$$ is susceptible is $$S(\hat{\mathbf{x}},0) \exp[-\kappa_v(1-\Theta(\hat{\mathbf{x}},t))]$$. So the probability a random neighbour is susceptible is given by   \begin{align*} \Phi_S &= \int_V S(\hat{\mathbf{x}},0)\int_\kappa \frac{\kappa P(\kappa)}{{\left\langle {\kappa} \right \rangle}} f(|\hat{\mathbf{x}}-\mathbf{x}|) e^{-\kappa(1-\Theta(\hat{\mathbf{x}},t))} \, \mathrm{d}\kappa \, \mathrm{d}\hat{\mathbf{x}}\\ &= \frac{\int_V S(\hat{\mathbf{x}},0) \Psi'(\Theta(\hat{\mathbf{x}},t)) f(|\hat{\mathbf{x}}-\mathbf{x}|) \, \mathrm{d}\hat{\mathbf{x}}}{{\left\langle {\kappa} \right \rangle}}. \end{align*} From this, $$\frac{\partial {}}{\partial {t}}{\Theta} = -\beta\Phi_I=-\beta(\Theta-\Phi_R-\Phi_S)$$ becomes   \begin{align*} \frac{\partial {}}{\partial {t}}{\Theta}(\mathbf{x},t) &= - \beta \Theta(\mathbf{x},t) + \gamma(1-\Theta(\mathbf{x},t))\\ &\quad + \beta \frac{\int_{V}S(\hat{\mathbf{x}},0) \Psi'(\Theta(\hat{\mathbf{x}},t)) f(|\hat{\mathbf{x}}-\mathbf{x}|) \, \mathrm{d}\hat{\mathbf{x}}}{{\left\langle {\kappa} \right \rangle}}. \end{align*} This derivation yields the following set of coupled equations for SIR propagation on RSNs:   \begin{align} \frac{\partial {}}{\partial {t}}{\Theta}(\mathbf{x},t) &= - \beta \Theta(\mathbf{x},t) + \gamma(1-\Theta(\mathbf{x},t)) \nonumber \label{eq:dtheta}\\ &\quad+ \beta \frac{\int_{V}S(\hat{\mathbf{x}},0)\Psi'(\Theta(\hat{\mathbf{x}},t)) f(|\hat{\mathbf{x}}-\mathbf{x}|) \, \mathrm{d}\hat{\mathbf{x}}}{{\left\langle {\kappa} \right \rangle}}\\ \end{align} (2a)  \begin{align} S(\mathbf{x},t) &= S(\mathbf{x},0)\Psi(\Theta(\mathbf{x},t)),\\ \end{align} (2b)  \begin{align} \frac{\partial {}}{\partial {t}}{R}(\mathbf{x},t) &= \gamma (1-S(\mathbf{x},t)-R(\mathbf{x},t))\\ \end{align} (2c)  \begin{align} I(\mathbf{x},t) &=1-S(\mathbf{x},t)-R(\mathbf{x},t) \label{eq:I} \end{align} (2d) with $$\Psi(\Theta) = \int_0^\infty e^{-\kappa(1-\Theta)} P(\kappa) \, \mathrm{d}\kappa$$. Our initial conditions are $$\Theta(\mathbf{x},0)=1$$ and $$R(\mathbf{x},0)=0$$, with $$S(\mathbf{x},0) = 1- I(\mathbf{x},0)$$. The equation for $$\Theta$$ is a nonlocal evolution equation with an integral over space and a derivative with respect to time. The nonlocal effects are captured through the convolution integral. To emphasize that the nonlocal interactions are captured by an integral, we refer to this as a partial integro-differential equation (PIDE): Unlike [35] we obtain an evolution equation for $${\Theta}(\mathbf{x},t)$$ that captures the RSN’s spatial structure in the integral term. We expect the equation for $$\Theta$$ to be similar to the Fisher–KPP equation for a spreading population [68, 69] $$\dot{u} = ru(1-u/K) + D u_{xx}$$, with the spatial integral playing the role of the nonlinear and diffusion terms. The spatial integral captures the network’s structure by including nonlocal interactions through $$f$$ and the degree distribution through $$\Psi$$. In Appendix B, we derive the Fisher–KPP equation as an approximation of (2a). Other nonlocal versions of the Fisher–KPP equation have been studied [70, 71], including some for disease spread [9, 39–41]. These are based on a mass-action mixing assumption and involve a convolution of a distance kernel with $$u$$ directly rather than a nonlinear function of $$u$$. Since stochastic simulations on networks are often hard to analyse and understand, this explicit analytic equation provides a powerful tool to study SIR dynamics on random networks with versatile spatial structure, fully taking into account the distribution of expected degrees and the distance kernel. In particular, the stochastic simulations to be discussed in Section 3 reveal the existence of travelling waves. In the next subsection, we use the PIDE to analytically derive properties of travelling waves in the analytic model, which, as we demonstrate in Section 3, accurately describe the travelling waves observed in the stochastic simulations. 2.3 Analytic properties of travelling waves on RSNs The PIDE formulation provides a powerful technical tool to derive precise quantitative insight in the spread of SIR disease on RSNs. As a major application, we can derive properties of travelling waves allowed by the PIDE. In particular, we derive an explicit expression for the wave speed of a 1D travelling wave, and identify conditions on the spatial decay of the kernel for the travelling wave to exist. This wave has a ‘pulled front’: That is, its speed is set by the nodes in the leading edge. We derive the wave speed assuming that the domain $$V$$ is the real line. We can write $$\Theta(x,t) = \theta(\xi(x,t))$$ where $$\xi(x,t) = x-ct$$ and $$c$$ is the wave speed. We will assume it is travelling from left to right. The wave travels into a region where $$S(x,t)\approx 1$$, $$I(x,t) \ll 1$$ and $$R(x,t) \ll 1$$. Very little transmission has occurred, so ahead of the travelling wave $$\Theta \approx 1$$. We assume $$\xi$$ is in the leading edge and write $$\Theta(x,t) = \theta(\xi(x,t))=1-\epsilon h(\xi(x,t))$$. In this leading edge of the wave $$\epsilon h(\xi) \ll 1$$, while in the bulk of the wave, $$\epsilon h(\xi)$$ may be comparable to $$1$$. We focus on the behaviour in the leading edge of the wave and transform equation (2a) for $$\frac{\partial {}}{\partial {t}}\Theta$$ into an equation for $$h$$ by expanding $$\Psi'(\theta(\xi))$$ as a Taylor Series about $$\theta=1$$:   \[ \Psi'(\theta) = \Psi'(1) - \epsilon h(\xi) \Psi''(1) + \frac{\epsilon^2 h(\xi)^2}{2}\Psi'''(1) + \cdots \] Note that $$\Psi^{(n)}(1) = {\left\langle {\kappa^n} \right \rangle}$$. Substituting this into the integral (and taking $$S(x,0)=1$$ ahead of the wave), we arrive at   \begin{align*} \epsilon c h'(\xi) &= -\beta(1-\epsilon h(\xi)) + \epsilon \gamma h(\xi) \\ &\quad +\beta \frac{\int_{-\infty}^\infty [\Psi'(1) - \epsilon h(\hat{\xi}) \Psi''(1)]\,f(|\hat{\xi}-\xi|) \, \mathrm{d}\hat{\xi}}{{\left\langle {\kappa} \right \rangle}}\\ & \quad + \beta \frac{\int_{-\infty}^\infty{\mathcal{O}}(\epsilon^2h(\hat{\xi})^2)\Psi'''(1)\,f(|\hat{\xi}-\xi|) \mathrm{d}\hat{\xi}}{{\left\langle {\kappa} \right \rangle}} \end{align*} Because $$\Psi'(1)={\left\langle {\kappa} \right \rangle}$$, the $${\mathcal{O}}(1)$$ terms cancel. We neglect the $${\mathcal{O}}(\epsilon^2 h(\hat{\xi})^2)$$ terms by arguing that if $$|\hat{\xi}-\xi|$$ is not large then $$\epsilon^2h(\hat{\xi})^2$$ is small (at the leading edge), and if $$|\hat{\xi}-\xi|$$ is large then $$f(|\hat{\xi}-\xi|)$$ is small (away from the leading edge). This leaves the linear homogeneous equation for $$h$$  \[ ch'(\xi) = (\beta+\gamma)h(\xi) - \beta \frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}\int_{-\infty}^\infty h(\hat{\xi})f(|\hat{\xi}-\xi|) \, \mathrm{d}\hat{\xi} \] We anticipate $$h(\xi) \approx e^{-\alpha \, \xi}$$ for some $$\alpha$$, yielding:   \[ -c\alpha = (\beta+\gamma) - \beta \frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}\int_{-\infty}^\infty e^{-\alpha (\hat{\xi}-\xi)} f(|\hat{\xi}-\xi|) \, \mathrm{d}\hat{\xi} \] Setting $$u = \hat{\xi}-\xi$$ the integral becomes $$\int_{-\infty}^\infty e^{-\alpha u} f(|u|)\, \mathrm{d}u$$. This is the bilateral Laplace transform of $$f(|x|)$$, which we denote $${\mathcal{L}}{}[\,f](\alpha)$$. Performing some further algebra yields   \begin{equation} \frac{c}{\beta+\gamma} = -\frac{1}{\alpha} + {\mathcal{R}_0} \frac{{\mathcal{L}}{}[\,f](\alpha)}{\alpha} \label{eq:alpha_c} \end{equation} (3) where $${\mathcal{R}_0} = \beta{\left\langle {\kappa^2} \right \rangle}/(\beta+\gamma){\left\langle {\kappa} \right \rangle}$$ (this is the basic reproductive number of the SIR disease on the network, which is the typical number of infections caused by an infected individual early in an epidemic [31]). There are infinitely many solution pairs $$\alpha,c$$. Following results for the Fisher–KPP equation we expect that the observed solution has the minimum wave speed. Setting $$\frac{\partial {c}}{\partial {\alpha}}=0$$ and performing some algebra shows that this occurs when   \begin{equation} \alpha{\mathcal{L}}{}[xf(|x|)](\alpha) + {\mathcal{L}}{}[\,f](\alpha) = \frac{1}{{\mathcal{R}_0}} \label{eqn:find_alpha} \end{equation} (4) We can solve this equation to find $$\alpha$$. Finally, substituting (3) into (4) gives   \begin{equation} c = -\beta\frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}} {\mathcal{L}}{}[xf(|x|)](\alpha) = -(\beta+\gamma) {\mathcal{R}_0} {\mathcal{L}}{}[xf(|x|)](\alpha) \label{eqn:find_c} \end{equation} (5) Note that the the wave speed $$c$$ and the spatial decay rate $$\alpha$$ depend on the distribution of expected degrees only through $${\left\langle {\kappa} \right \rangle}$$ and $${\left\langle {\kappa^2} \right \rangle}$$, specifically through the ratio $${\left\langle {\kappa^2} \right \rangle}/{\left\langle {\kappa} \right \rangle}$$. Figure 4 investigates the dependence of the speed on the degree distribution. Specifically, for different $$f$$ proportional to $$e^{-|x|^{1+\epsilon}}$$, we investigate how the wavespeed depends on $${\mathcal{R}_0} = \beta {\left\langle {\kappa^2} \right \rangle}/{\left\langle {\kappa} \right \rangle}$$. We normalize time so that $$\beta+\gamma=1$$. As $${\mathcal{R}_0}$$ grows, eventually $$\alpha$$ changes very slowly, and $$c$$ is approximately proportional to $${\mathcal{R}_0}$$. So increasing $${\mathcal{R}_0}$$ by modifying the distribution of $$\kappa$$ increases the wavespeed, and the speed tends to be proportional to $${\left\langle {\kappa^2} \right \rangle}/{\left\langle {\kappa} \right \rangle}$$. Fig. 4. View largeDownload slide Dependence of predicted wave properties on degree distribution. The predicted values of $$\alpha$$ and $$c$$ as a function of $${\mathcal{R}_0}$$ for $$f(x) \propto e^{-|x|^{1+\epsilon}}$$. We measure how the wavespeed depends on $${\mathcal{R}_0} = \frac{\beta}{\beta+\gamma}\frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}$$ with time normalized so that $$\beta+\gamma=1$$. Fig. 4. View largeDownload slide Dependence of predicted wave properties on degree distribution. The predicted values of $$\alpha$$ and $$c$$ as a function of $${\mathcal{R}_0}$$ for $$f(x) \propto e^{-|x|^{1+\epsilon}}$$. We measure how the wavespeed depends on $${\mathcal{R}_0} = \frac{\beta}{\beta+\gamma}\frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}$$ with time normalized so that $$\beta+\gamma=1$$. Note that if we chose another normalization of time, such as setting $$\gamma=1$$ our results suggest that increasing $${\mathcal{R}_0}$$ by increasing $$\beta$$ would mean that the wavespeed would scale like $$(1+\beta)\beta$$ for larger $$\beta$$. Note also that at $${\mathcal{R}_0}=1$$, both $$\alpha$$ and $$c$$ are $$0$$. If $$f$$ does not decay at least exponentially fast, then the Laplace transforms above do not exist. Thus we infer that if $$f$$ does not decay at least exponentially fast, then there is no coherent travelling wave solution in the $$\rho\to\infty$$ limit, consistent with [41, 72]. In practice, for a finite population if the long tail is observed by the transmissions, we see hop-and-spread dynamics, while if the long tail is not sampled by the transmission we may still see travelling wave behaviours in a given realization. This result for the propagation of SIR disease on random spatial networks described by Eq. (1) precisely reproduces the seminal result of Mollison [41] that for a simple 1D model of disease propagation in continuous space, travelling wave solutions can exist only for kernels that decay exponentially or faster. This is no surprise, because although the model used by Mollison (and, in a similar fashion, in [9]), was formulated in a simplified setting involving a convolution of a distance kernel with $$u$$ directly rather than a nonlinear function of $$u$$, it is still similar to our new PIDE. It is a powerful result, though, that we can analytically demonstrate this same behaviour, originally derived for simplified continuous space models, on the discrete networks generated by the RSN formula Eq. (1), in the limit of large node density. 2.4 Epidemic final size A further prediction for stochastic SIR spread on RSNs that can be computed analytically from the PIDE is the final size of an epidemic in a large population. As $$t \to \infty$$ the system approaches an equilibrium. By setting $$\frac{\partial {}}{\partial {t}}\Theta = 0$$, we have   \[ \Theta(\mathbf{x},\infty) = \frac{\gamma}{\beta+\gamma} + \frac{\beta}{\beta+\gamma} \frac{S(\hat{\mathbf{x}},0)\int_V \Psi'(\Theta(\hat{\mathbf{x}},\infty)) f(|\hat{\mathbf{x}}-\mathbf{x}|) \, \mathrm{d}\hat{\mathbf{x}}}{{\left\langle {\kappa} \right \rangle}} \] If $$\mathbf{x}$$ is far from the point of introduction, or the introduction is so small that we can approximate $$S(\mathbf{x},0)$$ as $$1$$ everywhere, we can treat $$\Theta(\mathbf{x},\infty)$$ as spatially homogeneous. Then   \begin{align} \Theta &= \frac{\gamma}{\beta+\gamma} + \frac{\beta}{\beta+\gamma} \frac{\Psi'(\Theta)}{{\left\langle {\kappa} \right \rangle}} \, , \label{eqn:theta_fs}\\ \end{align} (6)  \begin{align} S &= \Psi(\Theta) \, . \end{align} (7) $$\Theta=1$$, $$S=1$$ is always a solution, but if there is a solution $$\Theta$$ between $$0$$ and $$1$$, it is the appropriate one for an epidemic. It exists precisely when $${\mathcal{R}_0}>1$$. This is the same relation as for a random Chung–Lu network without spatial structure [35]. Interestingly, this means the epidemic final size in RSNs is independent of the distance kernel and depends only on disease properties and the degree distribution. 3. Numerical results We now apply the analytic models from Section 2 to analyse and understand stochastic simulations of SIR disease spread on RSNs. We first consider a simple Gaussian spatial kernel and demonstrate that the nonlocal evolution equations quantitatively predict stochastic simulations. The theory of Section 2 predicts the occurrence of travelling waves for the Gaussian kernel, and we do indeed observe travelling waves in the stochastic simulations on RSNs that are closely matched by numerical solutions of travelling waves in the PIDE. We also verify the analytic properties that were derived for the travelling waves in Section 2. We then consider a specific subclass of RSNs on which we investigate the interplay between short-range and long-range spatial connections in the SIR dynamics. 3.1 travelling waves and correspondence between stochastic SIR simulation and PIDE model We first compare numerical simulations of PIDE System (2) with stochastic SIR simulations on RSNs with a simple Gaussian kernel to demonstrate that these equations closely describe the stochastic dynamics in the limit of large $$N$$. In the process, we will identify travelling waves and verify the analytic properties we derived for them. We consider stochastic simulations on one-dimensional (1D) and two-dimensional (2D) RSNs generated using an algorithm based on the linear-time algorithm of [73] for Chung–Lu networks (see Appendix A). We simulate epidemic spread starting from a localized initial condition using a stochastic simulation algorithm from [32]. Figure 5 shows travelling waves revealed by the stochastic simulations in 1D and 2D. These retain their shape as they propagate. Figure 5 also shows numerical PIDE solutions of System (2) for the 1D and 2D network problems, demonstrating good agreement between PIDE solution and stochastic simulations on the 1D and 2D RSNs with Gaussian spatial kernel. Fig. 5. View largeDownload slide Comparison of stochastic simulation of SIR disease dynamics with numerical solution of analytic PIDE, on 1D and 2D RSNs. In both panels, results of stochastic simulation on networks with $$N=10^6$$ nodes are given by black dots. (Left) 1D PIDE solution is given by solid blue line (spatial resolution $$\Delta x = 10^{-4}$$). (Right) 2D PIDE solution is given by mesh surface (spatial resolution $$\Delta x = 1/512$$). Networks in both panels are generated using kernel $$f(r) = \phi_{0,0.01}(r)$$, where $$\phi_{\mu,\sigma}(x)$$ is the probability density function for the normal random variable with mean $$\mu$$ and standard deviation $$\sigma$$. All nodes have expected degree $$\kappa=20$$. Initial conditions: (left) 10% of nodes in the interval $$[0,0.01)$$ are initially infected, and (right) 10% of nodes in the square $$[0,1/32)\times[0,1/32)$$ are initially infected. Fig. 5. View largeDownload slide Comparison of stochastic simulation of SIR disease dynamics with numerical solution of analytic PIDE, on 1D and 2D RSNs. In both panels, results of stochastic simulation on networks with $$N=10^6$$ nodes are given by black dots. (Left) 1D PIDE solution is given by solid blue line (spatial resolution $$\Delta x = 10^{-4}$$). (Right) 2D PIDE solution is given by mesh surface (spatial resolution $$\Delta x = 1/512$$). Networks in both panels are generated using kernel $$f(r) = \phi_{0,0.01}(r)$$, where $$\phi_{\mu,\sigma}(x)$$ is the probability density function for the normal random variable with mean $$\mu$$ and standard deviation $$\sigma$$. All nodes have expected degree $$\kappa=20$$. Initial conditions: (left) 10% of nodes in the interval $$[0,0.01)$$ are initially infected, and (right) 10% of nodes in the square $$[0,1/32)\times[0,1/32)$$ are initially infected. Figure 6 compares the analytically predicted wave speed $$c$$ with wave speeds obtained from stochastic simulations. The small black circles of Fig. 6 give the wave speed observed in stochastic simulations. They converge to the analytic prediction (cyan dashed), but the convergence is very slow. Due to finite density, the exponentially decaying front is eventually truncated at its leading edge as shown in Fig. 7. This truncation slows the wave because the foremost infected nodes play a large role in setting the wave speed. This has been analysed for similar fronts in other stochastic systems, for which the leading order error in the wave speed decays proportionally to $$1/[\ln \rho]^2$$ [74]. More nodes are needed to improve the fit, as seen in Fig. 6. Figure 7 confirms that the stochastic simulation front and the numerically calculated PIDE front are nearly exponential with slope close to the predicted $$-\alpha$$. For the stochastic simulation (green dots), local densities in the exponentially decaying profile smaller than about 10$$^{-4}$$ cannot be represented because there are insufficient points in the simulation (the local densities effectively drop down to zero to the right of $$x\sim 0.5$$, and these zero values end up outside the range of the vertical logarithmic axis of the figure). By having larger $$\rho$$, more of the leading edge is observed resulting in wave speeds closer to the analytic prediction (Fig. 6). Figure 7 also shows that, in the numerical PIDE simulations, the exponential profile is truncated when the density of infected individuals approaches machine accuracy. Fig. 6. View largeDownload slide Comparison of wave speed for SIR disease dynamics observed in stochastic simulations (small black circles) and analytic prediction (dashed cyan line). For the stochastic simulations, we generate $$n_{rep} = 25$$ 1D spatial networks for each node density $$\rho$$ using a distance kernel $$f(|x|) = \phi_{0,0.01}(|x|)$$, where $$\phi_{\mu,\sigma}(x)$$ is the probability density function for the normal random variable with mean $$\mu$$ and standard deviation $$\sigma$$. All nodes have expected degree $$\kappa=10$$. For each network we realize one SIR simulation with disease parameters $$\beta=1$$ and $$\gamma=3$$. The black circles show the average wave speed resulting from 25 network realizations (vertical bars indicate 95% confidence intervals). The dotted black curve represents the expected convergence behaviour of $$c^*- K/[\ln \rho]^2$$ [74] to the analytically predicted wave speed $$c^*$$, where the constant $$K$$ is obtained by fitting the curve to the rightmost black circle, which has the smallest error bar since it corresponds to the largest node density $$\rho$$. Fig. 6. View largeDownload slide Comparison of wave speed for SIR disease dynamics observed in stochastic simulations (small black circles) and analytic prediction (dashed cyan line). For the stochastic simulations, we generate $$n_{rep} = 25$$ 1D spatial networks for each node density $$\rho$$ using a distance kernel $$f(|x|) = \phi_{0,0.01}(|x|)$$, where $$\phi_{\mu,\sigma}(x)$$ is the probability density function for the normal random variable with mean $$\mu$$ and standard deviation $$\sigma$$. All nodes have expected degree $$\kappa=10$$. For each network we realize one SIR simulation with disease parameters $$\beta=1$$ and $$\gamma=3$$. The black circles show the average wave speed resulting from 25 network realizations (vertical bars indicate 95% confidence intervals). The dotted black curve represents the expected convergence behaviour of $$c^*- K/[\ln \rho]^2$$ [74] to the analytically predicted wave speed $$c^*$$, where the constant $$K$$ is obtained by fitting the curve to the rightmost black circle, which has the smallest error bar since it corresponds to the largest node density $$\rho$$. Fig. 7. View largeDownload slide A comparison of stochastic simulation and numerical PIDE solution on a logarithmic scale, using the 1D solutions from Fig. 5. Both stochastic and numerical solutions experience leading edge truncation. For the stochastic simulation it is due to the finite number of nodes, while for the numerical solution it is due to numerical precision ($$\sim 10^{-16}$$). The slope of the leading edge is close to the asymptotic prediction of $$-\alpha$$. Fig. 7. View largeDownload slide A comparison of stochastic simulation and numerical PIDE solution on a logarithmic scale, using the 1D solutions from Fig. 5. Both stochastic and numerical solutions experience leading edge truncation. For the stochastic simulation it is due to the finite number of nodes, while for the numerical solution it is due to numerical precision ($$\sim 10^{-16}$$). The slope of the leading edge is close to the asymptotic prediction of $$-\alpha$$. 3.2 A subclass of RSNs with long-range connections We now consider a specific subclass of RSNs with short-range and long-range spatial connections on which we will investigate the effect long-range connections have on SIR dynamics. The networks in this subclass have an average degree of $$k$$, with edges distributed such that nodes within some distance $$r_0$$ are very likely to be connected, and nodes of greater distance are very unlikely to be connected. These are modelled after the Watts–Strogatz small-world networks [47] with the number of local edges decreasing as the long-range connections increase so that the average degree remains fixed. Specifically, we set $$\kappa = k$$ for all nodes, choose $$r_0$$, and define a reference node density $$N_0:=k/(\pi r_0^2)$$. We place $$N \ge N_0$$ nodes uniformly at random into the unit torus $$[0,1]\times[0,1]$$ with periodic boundaries and $$|V|=1$$. We assume $$r_0<1/2$$ so that disks of radius $$r_0$$ easily fit within the torus. In particular this forces $$\pi r_0^2<1$$. We choose a small number $$\epsilon$$ and define the spatial kernel as   \begin{equation} f(d_{uv}) = \begin{cases} \frac{N_0}{k}\left[1-\epsilon\frac{1-\pi r_0^2}{\pi r_0^2}\right]& d_{uv} < r_0\\ \frac{N_0}{k}\epsilon & d_{uv} \geq r_0 \end{cases}. \label{eq:kernel} \end{equation} (8) Note that $$\epsilon$$ cannot be greater than $$\pi r_0^2/(1-\pi r_0^2)$$. Equation (1) specializes to   \begin{equation} p_{uv} = \min\left(k \frac{f(d_{uv})}{N}, 1\right) = \frac{k f(d_{uv})}{N}, \label{eq:pspec} \end{equation} (9) since, for small $$\epsilon$$, $$f(d_{uv}) \le N_0/k$$ so $$kf(d_{uv})/N \le 1$$, assuming $$\pi r_0^2<1$$. Then   \begin{equation} p_{uv} = \begin{cases} \frac{N_0}{N} \left(1- \epsilon\frac{1-\pi r_0^2}{\pi r_0^2}\right) & d_{uv} <r_0\\ \frac{N_0}{N} \epsilon & d_{uv}\geq r_0 \end{cases}. \label{eq:puv} \end{equation} (10) The expected number of nodes within distance $$r_0$$ is $$\pi r_0^2 N$$ and beyond distance $$r_0$$ is $$(1-\pi r_0^2)N$$. The total number of expected neighbours in the graph is $$\pi r_0^2 N \frac{N_0}{N} \left(1-\epsilon(1-\pi r_0^2)/(\pi r_0^2)\right)+(1-\pi r_0^2)N \frac{N_0}{N} \epsilon=k$$. When we place $$N=N_0$$ nodes and set $$\epsilon=0$$, all nodes within a distance $$r_0$$ are connected ($$p_{uv}=1$$) and no other nodes are joined. This is a random geometric graph [75]. When $$\epsilon=0$$ and $$N$$ is chosen greater than $$N_0$$, there are more than $$k$$ nodes within a distance $$r_0$$, but not all of these nodes are connected ($$p_{uv} < 1$$), such that the expected degree is still $$k$$. For any choice of $$N \geq N_0$$, as $$\epsilon$$ increases, the number of long-range connections grows proportionally to $$\epsilon$$ with a corresponding decrease in the short-range connections. In what follows we consider networks of this subclass with $$N>N_0$$ and varying $$\epsilon$$, i.e., a varying fraction of long-range connections. 3.3 An example of SIR disease propagation on RSNs with short-range and long-range connections: travelling waves and hop-and-spread dynamics We now perform stochastic simulations [32] of SIR disease on RSNs with spatial kernel (8) and $$k=20$$, for three values of $$\epsilon$$. We choose $$N_0= 10^5$$ such that $$r_0=\sqrt{k/(\pi N_0)} = 8.0\times10^{-3}$$, and place $$N=10^6$$ random nodes. Nodes transmit infection with rate $$\beta=1$$ and recover with rate $$\gamma=3$$. The simulated behaviour shown in Fig. 8 reveals: When there are few long-range connections ($$\epsilon \to 0$$), the disease spreads in a wave-like, coherent spatial pattern outwards from the source. With an increased number of long-range connections, we see hop-and-spread dynamics. The spread is locally wave-like until a long-range transmission occurs. It then generates a new focal region. Once the number of infections is large enough that the disease experiences its first long-range transmission, it usually starts causing many long-range transmissions. As the number of long-range transmissions increases, the disease spreads uniformly throughout the population. Sufficient long-range connections mean the epidemic behaves like it would in a globally-mixed population. Fig. 8. View largeDownload slide Stochastic simulation of SIR disease dynamics on a 2D RSN with spatial kernel (8) with $$N=10^6$$ nodes. The fraction of infected nodes, $$I$$, is shown, with a colour scale that ranges from blue ($$I=0$$) to yellow ($$I=0.5$$). The density of infected nodes at time $$t=10$$ seconds is shown for networks with (left) highly local spatial structure ($$\epsilon=10^{-10}$$), (middle) intermediate spatial structure ($$\epsilon=10^{-8.25}$$) and (right) no local spatial structure ($$\epsilon=\pi r_0^2$$, i.e., $$f(d_{uv})=1$$, uniform spatial kernel). Initial condition: the 1000 nodes closest to (0.5,0.5) are initially infected. The waves in the first two panels behave as travelling waves, as can be seen in the supplementary movies 1–3 (one for each panel). Fig. 8. View largeDownload slide Stochastic simulation of SIR disease dynamics on a 2D RSN with spatial kernel (8) with $$N=10^6$$ nodes. The fraction of infected nodes, $$I$$, is shown, with a colour scale that ranges from blue ($$I=0$$) to yellow ($$I=0.5$$). The density of infected nodes at time $$t=10$$ seconds is shown for networks with (left) highly local spatial structure ($$\epsilon=10^{-10}$$), (middle) intermediate spatial structure ($$\epsilon=10^{-8.25}$$) and (right) no local spatial structure ($$\epsilon=\pi r_0^2$$, i.e., $$f(d_{uv})=1$$, uniform spatial kernel). Initial condition: the 1000 nodes closest to (0.5,0.5) are initially infected. The waves in the first two panels behave as travelling waves, as can be seen in the supplementary movies 1–3 (one for each panel). The coherent spatial spread and globally-mixed regimes behave largely deterministically, but the intermediate hop-and-spread regime is dominated by stochastic effects. In a large enough domain (or equivalently, for small enough $$r_0$$), an epidemic spreading with any $$\epsilon>0$$ will eventually exhibit hop-and-spread dynamics. How long the initial wave travels before the hop-and-spread dynamics begin depends on the time until the disease crosses a long-range connection. We can gain insight by initially ignoring the long-range connections. Then the disease spreads in an annulus. The width and rate of spread of the annulus is roughly constant as the radius increases. Thus the number of infected nodes increases linearly. From the frequency of long-range edges and the transmission probability per edge, we can estimate the expected number of infections until a long-range transmission occurs. Once one long-range transmission occurs, the number infected is likely high enough that many other long-range transmissions begin to occur. So the typical dynamics have an early phase in which local transmissions dominate, yielding a growing localized annulus of infection. Then long-range transmissions followed by new foci occur, and this dominates the later dynamics. 3.4 Factors affecting accuracy of analytic approximation The analytic equations are effectively a continuum approximation to the discrete dynamics of the process. To be a valid approximation, it needs to be sensible to talk about the density of nodes and edges near a point in space. That is, we need the transition from Question 1 to Question 2 in section 2.2 to be appropriate. Thus we need a separation of scales: we need a scale which is small compared with the dynamics of the spatial spread, but large enough that the local expected behaviour gives a very accurate prediction of the number of nodes within the region having a given status. The specific conditions under which the analytic equations become appropriate are both disease and network-dependent. Let us consider a region $$X$$ which is small compared with the spatial scales over which the disease spread is predicted to vary. The predictions will tend to be inaccurate when they predict a very small number of infected individuals in $$X$$. We can classify these into four situations If the local reproductive number is less than one, the inaccuracies are not important because the disease is dying off in both the prediction and the realization. If $$X$$ is slightly ahead of an advancing front, the error will tend to cause a reduction in the front’s speed, but otherwise the predictions perform well. If $$X$$ is far ahead of any front, then the impact may be significant: an introduced single case where the expected number of cases is much less than $$1$$ can lead to a new focus of spread, as seen in Fig. 8. If the density of nodes is low enough that $$X$$ does not have many nodes or clustering is significant, then the arguments leading to the governing equations break down. The first two issues lead to relatively small errors, while the final two can lead to significant errors. 4. Discussion and conclusion In this article, we have derived new analytic models of PIDE type for the spread of SIR disease on random networks with spatial structure. We derive these models for a class of RSNs that intrinsically incorporates spatial structure in a way that crucially remains analytically tractable. These models are important to understand and quantify the dynamics of SIR disease spread on networks with spatial structure. We observe nonlinear travelling waves on RSNs and are able to precisely characterize these waves analytically. This is the first quantitatively accurate analytic description of nonlinear travelling waves on random graphs with spatial structure. Our SIR disease model on spatial graphs extends analytic approaches that were previously only available to describe dynamics at the level of populations or on networks without spatial structure, to the real-life case of networks characterized by important spatial structure. This provides a theoretical modelling framework, for example, for recent observations of how epidemics like Ebola evolve in modern connected societies, with long-range connections seeding new focal points from which the epidemic spreads locally (see [11], in particular Video 1). Figure 9 shows the observed hops seen from virus genome sequencing for the 2013-2016 Ebola outbreak [11], which can be built into RSN models. Fig. 9. View largeDownload slide Observed hop distances inferred from viral genome sequencing for the 2013–2016 Ebola outbreak in West-Africa. Data courtesy of Rambaut [11]. The estimated probability density is obtained from raw binned data by Kernel Density Estimation using a Gaussian kernel function. The dotted lines denote the 50th and 95th percentiles. Fig. 9. View largeDownload slide Observed hop distances inferred from viral genome sequencing for the 2013–2016 Ebola outbreak in West-Africa. Data courtesy of Rambaut [11]. The estimated probability density is obtained from raw binned data by Kernel Density Estimation using a Gaussian kernel function. The dotted lines denote the 50th and 95th percentiles. There are many possibilities for further applications in disease modelling. For example, the SIR travelling wave models can be applied to realistic simulations of historic epidemics such as the plague in medieval Europe [9], incorporating in a precise manner the geography and estimated historical population density maps. Further potential applications include models for a diversity of areas such as neuronal networks in the human brain, spread of animal and plant species in fragmented habitats, and propagation of civil unrest, or public health epidemics such as obesity and smoking, in human societies. The RSN class is also amenable to graph theoretic analysis. For example, it would be of great interest to investigate thresholds for the existence of a giant component for RSNs with various distance kernels. This was successfully done for non-spatial random graphs with given degree sequences by Reed et al. [36, 76]. The techniques in [77] show that this should be straightforward for RSNs. Also, the ODE systems resulting from applying the analytic edge-based compartmental modelling to SIR disease spread on networks without spatial structure [35] have been shown in a mathematically rigorous way to be exact approximations of the stochastic SIR dynamics in the limit of large node size. Our analytic PIDE models for SIR spread on networks with spatial structure were derived using similar assumptions, and the close match we obtain between stochastic and PIDE simulations as, e.g., in Fig. 5, indicates that there is a mathematical conjecture that our PIDE models are also exact in the limit of large node density. It may be possible to prove this mathematically using techniques similar to the ones used in [37, 38]. Furthermore, the RSN model defined by (1) is versatile as it allows flexible degree distributions and distance kernels. There are some further generalizations that were not considered in this article but offer compelling prospects for building realistic networks models that remain analysable. An obvious generalization is to allow $$\rho(\mathbf{x})$$ to describe an inhomogeneous spatial density of nodes. This could, for example, be used to model different population densities in cities and rural areas in the context of realistic spatial disease spreading models, or different densities of neurons in different parts of the brain. Similarly, instead of letting the distance kernel depend on nodal distances $$d_{uv}$$, one can consider more general kernels $$f(\mathbf{x}_u,\mathbf{x}_v)$$ that directly depend on the coordinates of the nodes. For example, this could model people living in cities preferentially connecting to people in the same and other cities, while connections within and between rural individuals could follow different patterns. This type of modelling is especially compelling in the era of big data where real data may be used to build analysable spatial random network models that faithfully mirror real-world spatial networks. For example, [10] studies the interplay between short-scale commuting flows and long-range airline traffic in shaping the spatiotemporal pattern of a global epidemic, and [11] found that during the 2013–2016 Ebola outbreak viral lineages moved according to a classic ‘gravity’ model, with more intense migration between larger and more proximate population centres. Using empirical data or inferred spatial connectivities, all these kinds of effects can be incorporated in our RSN models of (1), with clear potential for highly realistic models of disease propagation on spatial networks that retain the analytic tractability of the approach. Funding Joel C. Miller was funded by the Global Good Fund through the Institute for Disease Modeling and by a Larkins Fellowship from Monash University. Jamieson L. Kaiser was funded by an Australian Mathematical Sciences Institute Vacation Research Scholarship. John C. Lang acknowledges funding from the Natural Sciences and Engineering Research Council of Canada. Appendix A. Fast network generation At first sight, generating networks from the RSN class defined by Eq. (1) appears to be an $${\mathcal{O}}(N^2)$$ operation as each of the $$\binom{N}{2}$$ edges exists independently of the others. However, by modifying methods developed to generate Erdős–Rényi and Chung–Lu networks [44, 78] in linear time, we arrive at a much more efficient process. This makes large RSNs practical for simulation and analysis. We briefly outline the method, under the assumption that the distance kernel $$f$$ is bounded above and is monotonically decreasing. We divide the space $$V$$ into regions and order the nodes in each region by decreasing expected degree. We consider a region $$X_i$$, and define $$u$$ to be the first node in that region. The probability that $$u$$ will have an edge with any subsequent node is bounded above by $$q_0=\min(1,\kappa_u^2 f_{\text{max}}/{\left\langle {\kappa} \right \rangle})$$ where $$f_{\text{max}}= f(0)$$ is the maximum of $$f$$. We then choose a number $$r$$ from a geometric distribution with success probability $$q_0$$. So $$r$$ is chosen from a negative binomial distribution giving the number of unsuccessful trials (each having success probability $$q_0$$) before the first success. We set $$v_1$$ to be the $$r$$th node following $$u$$ in $$X_i$$. This is equivalent to $$v_1$$ being the first node chosen when sequentially considering each node with probability $$q_0$$. Thus the ‘candidate neighbour’ $$v_1$$ was chosen with probability $$q_0$$ which is at least $$p_{uv_1}=\min(1,\kappa_u \kappa_{v_1} f(d_{uv})/{\left\langle {\kappa} \right \rangle})$$. It is possible that $$v_1$$ is a “false positive”. To decide this, we generate a new random number uniformly from $$(0,1)$$, and if it is less than $$p_{uv_1}/q_0$$, we decide that $$v_1$$ was correctly chosen, and we add an edge between them. Otherwise $$v_1$$ is a false positive and we do not add the edge. We then enter an iterative step. After processing $$v_i$$, we set $$q_i = \min (1,\kappa_u \kappa_{v_i} f_{\text{max}}/{\left\langle {\kappa} \right \rangle})$$ and jump to the next candidate neighbour $$v_{i+1}$$. Again $$v_{i+1}$$ may be a false positive, and we place an edge between $$u$$ and $$v_{i+1}$$ with probability $$q_{uv_{i+1}}/q_i$$ where $$q_{uv_{i+1}} = \min (1,\kappa_u \kappa_{v_{i+1}} f(d_{uv_{i+1}})/{\left\langle {\kappa} \right \rangle})$$. As the iteration progresses, $$\kappa_{v_i}$$ decreases so $$q_i$$ decreases, and the jumps between successive $$v$$s become larger. The edges within each region are assigned in linear time. We next place edges between the node $$u$$ and other regions $$X_j$$. We find the minimum distance between $$u$$ and $$X_j$$ and use it to define $$f_{\text{max}}$$. We take $$w$$ to be the first node in $$X_j$$. We define $$q_0 = \min \{ 1, \kappa_u \kappa_w f_{\text{max}}/{\left\langle {\kappa} \right \rangle}\}$$. We choose $$v_1$$ from $$X_j$$ as before, and iterate. For $$X_j$$ farther from $$u$$ the jumps are larger. Appendix B. Relation to the Fisher–KPP equation We note that, under approximations that are appropriate in the case of a localized spatial kernel with $$\Theta$$ varying slowly in space we can convert (2) into the Fisher–KPP equation $$u_t = ru(1-u/K) + D u_{xx}$$. We demonstrate this in the 1D case. We start by assuming that $$S(x,0)$$ is approximately $$1$$. We then assume that $$f$$ is sufficiently localized and $$\Theta$$ varies sufficiently slowly that we can expand $$\Psi'(\Theta(x,t))$$ as a Taylor Series in $$x$$ to third order.   \begin{align*} \int_{-\infty}^\infty \Psi'(\Theta(\hat{x},t)) f(|\hat{x}-x|) \, \mathrm{d}\hat{x} = & \int_{-\infty}^\infty \Psi'(\Theta(x,t)) \, \mathrm{d}\hat{x}\\ & + \int_{-\infty}^\infty (\hat{x}-x)\frac{\partial {}}{\partial {x}}\Psi'(\Theta(x,t)) f(|\hat{x}-x|)\, \mathrm{d}\hat{x} \\ &+ \int_{-\infty}^\infty \frac{(\hat{x}-x)^2}{2}\left(\frac{\partial^{2} {}}{\partial x^{2}} \Psi'(\Theta(x,t)\right) f(|\hat{x}-x|) \, \mathrm{d}\hat{x} \\ &+ \cdots\\ \approx & \Psi'(\Theta(x,t)) + C\frac{\partial^{2} {}}{\partial x^{2}}\Psi'(\Theta(x,t)) \end{align*} where $$C = \int_{-\infty}^\infty y^2f(|y|)/2 \, \mathrm{d}y$$, and we have used the fact that $$\int_{-\infty}^\infty y f(|y|) \, \mathrm{d}y=0$$ by symmetry. So taking $$S(x,0)=1$$, our equation for $$\dot{\Theta}$$ is   \[ \dot{\Theta} = - \beta \Theta + \gamma(1-\Theta) + \beta \frac{\Psi'(\Theta)}{{\left\langle {\kappa} \right \rangle}} + \frac{\beta C}{{\left\langle {\kappa} \right \rangle}}\frac{\partial^{2} {}}{\partial x^{2}} \Psi'(\Theta) \] Now setting $$u = 1-\Theta$$ and assuming $$u$$ is small, we use further Taylor expansions and the fact that $$\Psi^{(n)}(1) = {\left\langle {\kappa^n} \right \rangle}$$ to find   \begin{align*} \dot{u} &= \beta (1-u) - \gamma u - \beta \frac{\Psi'(1-u)}{{\left\langle {\kappa} \right \rangle}} - \frac{\beta C}{{\left\langle {\kappa} \right \rangle}}\frac{\partial^{2} {}}{\partial x^{2}} \Psi'(1-u)\\ &\approx \beta (1-u) - \gamma u - \beta + \beta\frac{u{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}} - \beta \frac{u^2 {\left\langle {\kappa^3} \right \rangle}}{2{\left\langle {\kappa} \right \rangle}} - \frac{\beta C}{{\left\langle {\kappa} \right \rangle}}\frac{\partial^{2} {}}{\partial x^{2}}\Psi'(1-u)\\ &\approx \left(\beta \frac{{\left\langle {\kappa^2} \right \rangle}}{{\left\langle {\kappa} \right \rangle}}-\beta-\gamma \right) u - \beta \frac{{\left\langle {\kappa^3} \right \rangle}}{2{\left\langle {\kappa} \right \rangle}}u^2 - \frac{\beta C}{{\left\langle {\kappa} \right \rangle}}\frac{\partial^{2} {}}{\partial x^{2}} \left({\left\langle {\kappa} \right \rangle} - u {\left\langle {\kappa^2} \right \rangle}\right)\\ &= r u(1-u/K) + D\frac{\partial^{2} {}}{\partial x^{2}} u \end{align*} for appropriately chosen $$r$$, $$K$$, and $$D$$. So the Fisher–KPP equation arises out of an approximation of our PIDE on RSNs. In general these approximations are not accurate when, e.g., $$\Theta$$ is not close to $$1$$, the variation in $$\Theta$$ is fast enough that the Taylor Series expansions are poor approximations, or the spatial kernel is not localized. However, this relation does suggest that behaviours found for the Fisher–KPP equation are likely to occur for disease spread in our RSNs as well. Footnotes 1 Technically the number of neighbours is binomially distributed with parameters $$N-1$$ and $$p=\kappa_u/(N-1)$$ which can be treated as Poisson-distributed with mean $$\kappa_u$$ for large $$N$$, or if we take $$\rho$$ as fixed and place nodes as a true Poisson process, then $$N$$ is a random variable, and the number of neighbours truly is Poisson distributed. References 1. Grenfell B. T., Bjørnstad O. N. & Kappey J. ( 2001) Travelling waves and spatial hierarchies in measles epidemics. Nature , 414, 716– 723. Google Scholar CrossRef Search ADS PubMed  2. Childs J. E., Curns A. T., Dey M. E., Real L. A., Feinstein l., Bjørnstad O. N. & Krebs J. W. ( 2000) Predicting the local dynamics of epizootic rabies among raccoons in the United States. Proc. Natl. Acad. Sci. U.S.A , 97, 13666– 13671. Google Scholar CrossRef Search ADS PubMed  3. Biek R., Henderson J. C., Waller L. A., Rupprecht C. E. & Real L. A. ( 2007) A high-resolution genetic signature of demographic and spatial expansion in epizootic rabies virus. Proc. Natl. Acad. Sci. U.S.A , 104, 7993– 7998. Google Scholar CrossRef Search ADS PubMed  4. Steck F. & Wandeler A. ( 1980) The epidemiology of fox rabies in europe. Epidemiol. Rev. , 2, 71– 96. Google Scholar CrossRef Search ADS PubMed  5. Baer G. M. ( 1991) The Natural History of Rabies. 2nd edn, Chapter 13-titled as ’fox rabies’  ( Baer G. M. ed.). Boca Raton, Ann Arbor and Boston: CRC Press. 6. Benavides J. A., Valderrama W. & Streicker D. G. ( 2016) Spatial expansions and travelling waves of rabies in vampire bats. In Proc. R. Soc. B , vol. 283, UK: The Royal Society, pp. 20160328. Google Scholar CrossRef Search ADS   7. Bos K. I., Schuenemann V. J., Golding G. B., Burbano H. A., Waglechner N., Coombes B. K., McPhee J. B., DeWitte S. N., Meyer M., Schmedes S. et al.   ( 2011) A draft genome of yersinia pestis from victims of the black death. Nature , 478, 506– 510. Google Scholar CrossRef Search ADS PubMed  8. Christakos G., Olea R. A. & Yu H.-L. ( 2007) Recent results on the spatiotemporal modelling and comparative analysis of black death and bubonic plague epidemics. Public Health , 121, 700– 720. Google Scholar CrossRef Search ADS PubMed  9. Marvel S. A., Martin T., Doering C. R., Lusseau D. & Newman M. E. J. ( 2013) The small-world effect is a modern phenomenon. Preprint arXiv:1310.2636. 10. Balcan D., Colizza V., Gonalves B., Hu H., Ramasco J. J. & Vespignani A. ( 2009) Multiscale mobility networks and the spatial spreading of infectious diseases. Proc. Natl. Acad. Sci. U.S.A , 106, 21484– 21489. Google Scholar CrossRef Search ADS PubMed  11. Dudas G., Carvalho L. M., Bedford T., Tatem A. J., Baele G., Faria N. R., Park D. J., Ladner J. T., Arias A., Asogun D., et al.   ( 2017) Virus genomes reveal factors that spread and sustained the ebola epidemic. Nature , 544, 309– 315. Google Scholar CrossRef Search ADS PubMed  12. Fine P. E. M. & Carneiro I. A. M. ( 1999) Transmissibility and persistence of oral polio vaccine viruses: implications for the global poliomyelitis eradication initiative. Am. J. Epidemiol. , 150, 1001– 1021. Google Scholar CrossRef Search ADS PubMed  13. Kroiss S. J., Famulare M., Lyons H., McCarthy K. A., Mercer L. D. & Chabot-Couture G. ( 2017) Evaluating cessation of the type 2 oral polio vaccine by modeling pre-and post-cessation detection rates. Vaccine , 35, 5674– 5681. Google Scholar CrossRef Search ADS PubMed  14. Newman M. E. J. ( 2002) Spread of epidemic disease on networks. Phys. Rev. E , 66, 016128. Google Scholar CrossRef Search ADS   15. Meyers L. A. ( 2007) Contact network epidemiology: Bond percolation applied to infectious disease prediction and control. Bull. Am. Math. Soc. , 44, 63– 86. Google Scholar CrossRef Search ADS   16. Pastor-Satorras R., Castellano C., Van Mieghem P. & Vespignani A. ( 2015) Epidemic processes in complex networks. Rev. Mod. Phys. , 87, 925. Google Scholar CrossRef Search ADS   17. Bullmore E. & Sporns O. ( 2009) Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. , 10, 186– 198. Google Scholar CrossRef Search ADS PubMed  18. O’Dea R., Crofts J. J. & Kaiser M. ( 2013) Spreading dynamics on spatially constrained complex brain networks. J. R. Soc. Interface , 10, 20130016. Google Scholar CrossRef Search ADS PubMed  19. Ercsey-Ravasz M., Markov N. T., Lamy C., Van Essen D. C., Knoblauch K., Toroczkai Z. & Kennedy H. ( 2013) A predictive network model of cerebral cortical connectivity based on a distance rule. Neuron , 80, 184– 197. Google Scholar CrossRef Search ADS PubMed  20. Horvát S., Gămănu R., Ercsey-Ravasz M., Magrou L., Gămănu B., Van Essen D. C., Burkhalter A., Knoblauch K., Toroczkai Z. & Kennedy H. ( 2016) Spatial embedding and wiring cost constrain the functional layout of the cortical network of rodents and primates. PLoS Biol. , 14, e1002512. Google Scholar CrossRef Search ADS PubMed  21. Barthélemy M. ( 2011) Spatial networks. Phys. Rep. , 499, 1– 101. Google Scholar CrossRef Search ADS   22. Auffinger A., Damron M. & Hanson J. ( 2017) 50 years of first passage percolation. Providence, RI: American Mathematical Society. 23. Chung F. & Lu L. ( 2002) Connected components in random graphs with given expected degree sequences. Ann. Comb. , 6, 125– 145. Google Scholar CrossRef Search ADS   24. Newman M. E. J. ( 2003) The structure and function of complex networks. SIAM Rev. , 45, 167– 256. Google Scholar CrossRef Search ADS   25. Centola D., Eguíluz V. M. & Macy M. W. ( 2007) Cascade dynamics of complex propagation. Phys. A Stat. Mech. Appl. , 374, 449– 456. Google Scholar CrossRef Search ADS   26. Janssen H.-K., Müller M. & Stenull O. ( 2004) Generalized epidemic process and tricritical dynamic percolation. Phys. Rev. E , 70, 026114. Google Scholar CrossRef Search ADS   27. Watts D. J. ( 2002) A simple model of global cascades on random networks. Proc. Natl. Acad. Sci. U.S.A , 99, 5766– 5771. Google Scholar CrossRef Search ADS PubMed  28. Miller J. C. ( 2016) Complex contagions and hybrid phase transitions. J. Complex Netw. , cnv021. 29. Melnik S., Hackett A., Porter M. A., Mucha P. J. & Gleeson J. P. ( 2011) The unreasonable effectiveness of tree-based theory for networks with clustering. physical review E , 83, 036112. Google Scholar CrossRef Search ADS   30. Kermack W. O. & McKendrick A. G. ( 1927) A contribution to the mathematical theory of epidemics. In Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences , vol. 115. The Royal Society, pp. 700– 721. 31. Anderson R. M. & May R. M. ( 1991) Infectious Diseases of Humans . Oxford: Oxford University Press. 32. Kiss I. Z., Miller J. C. & Simon P. L. ( 2017) Mathematics of Epidemics on Networks: From Exact to Approximate Models . IAM. New York: Springer. Google Scholar CrossRef Search ADS   33. Abuelezam N. N., Rough K. & Seage III G. R. ( 2013) Individual-based simulation models of HIV transmission: Reporting quality and recommendations. PLoS One , 8, e75624. Google Scholar CrossRef Search ADS PubMed  34. Eaton J. W., Johnson L. F., Salomon J. A., Bärnighausen T., Eran B., Bershteyn A., Bloom D. E., Cambiano V., Fraser C., Hontelez J. A. C., et al.   ( 2012) HIV treatment as prevention: systematic comparison of mathematical models of the potential impact of antiretroviral therapy on HIV incidence South Africa. PLoS Med. , 9, e1001245. Google Scholar CrossRef Search ADS PubMed  35. Miller J. C., Slim A. C. & Volz E. M. ( 2012) Edge-based compartmental modelling for infectious disease spread. J. R. Soc. Interface , 9, 890– 906. Google Scholar CrossRef Search ADS PubMed  36. Molloy M. & Reed B. ( 1995) A critical point for random graphs with a given degree sequence. Random Struct. & Algorithms , 6, 161– 180. Google Scholar CrossRef Search ADS   37. Janson S., Luczak M. & Windridge P. ( 2014) Law of large numbers for the SIR epidemic on a random graph with given degrees. Random Struct. & Algorithms , 45, 724– 761. Google Scholar CrossRef Search ADS   38. Decreusefond L., Dhersin J.-S., Moyal P. & Tran V. C. ( 2012) Large graph limit for an SIR process in random network with heterogeneous connectivity. Ann. Appl. Prob. , 22, 541– 575. Google Scholar CrossRef Search ADS   39. Diekmann O. ( 1979) Run for your life. A note on the asymptotic speed of propagation of an epidemic. J. Differ. Equ. , 33, 58– 73. Google Scholar CrossRef Search ADS   40. Diekmann O. ( 1978) Thresholds and travelling waves for the geographical spread of infection. J. Math. Biol. , 6, 109– 130. Google Scholar CrossRef Search ADS PubMed  41. Mollison D. ( 1972) The rate of spatial propagation of simple epidemics. In Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Volume 3: Probability Theory . Berkeley and Los Angeles: The Regents of the University of California. 42. Bollobás B., Janson S. & Riordan O. ( 2007) The phase transition in inhomogeneous random graphs. Random Struct. & Algorithms , 31, 3– 122. Google Scholar CrossRef Search ADS   43. Norros I. & Reittu H. ( 2006) On a conditionally Poissonian graph process. Adv. Appl. Prob. , 38, 59– 75. Google Scholar CrossRef Search ADS   44. Miller J. C. & Hagberg A. ( 2011) Efficient generation of networks with given expected degrees. Proceedings of the 8th International Workshop on Algorithms and Models for the Web Graph . Berlin: Springer, pp. 115– 126. Google Scholar CrossRef Search ADS   45. Söderberg B. ( 2002) General formalism for inhomogeneous random graphs. Physical Review E , 66, 066121. Google Scholar CrossRef Search ADS   46. Gilbert E. N. ( 1959) Random graphs. Ann. Math. Stat. , 30, 1141– 1144. Google Scholar CrossRef Search ADS   47. Watts D. J. & Strogatz S. H. ( 1998) Collective dynamics of small-world networks. Nature , 393, 440– 442. Google Scholar CrossRef Search ADS PubMed  48. Newman M. E. J. & Watts D. J. ( 1999) Renormalization group analysis of the small-world network model. Phy. Lett. A , 263, 341– 346. Google Scholar CrossRef Search ADS   49. Kleinberg J. M. ( 2002) Small-world phenomena and the dynamics of information. In Advances in neural information processing systems ., vol. 15. Boston: MIT Press, pp. 431– 438. 50. Waxman B. M. ( 1988) Routing of multipoint connections. IEEE J. Selected Areas Commun. , 6, 1617– 1622. Google Scholar CrossRef Search ADS   51. Newman C. M. & Schulman L. S. ( 1986) One dimensional $$1/| j- i|^s$$ percolation models: The existence of a transition for $$s \leq 2$$. Commun. Math. Phys. , 104, 547– 571. Google Scholar CrossRef Search ADS   52. Biskup M. ( 2004) On the scaling of the chemical distance in long-range percolation models. Ann. Prob. , 32, 2938– 2977. Google Scholar CrossRef Search ADS   53. Trapman P. ( 2010) The growth of the infinite long-range percolation cluster. Ann. Prob. , 38, 1583– 1608. Google Scholar CrossRef Search ADS   54. Deijfen M., van der Hofstad R. & Hooghiemstra G. ( 2013) Scale-free percolation. In Annales de l’Institut Henri Poincaré, Probabilités et Statistiques , vol. 49. Institut Henri Poincaré, pp. 817– 838. Google Scholar CrossRef Search ADS   55. Deprez P. & Wüthrich M. V. ( 2015) Networks, random graphs and percolation. In Theoretical Aspects of Spatial-Temporal Modeling . Tokyo: Springer, pp. 95– 124. Google Scholar CrossRef Search ADS   56. Vladimirov N., Traub R. D. & Tu Y. ( 2011) Wave speed in excitable random networks with spatially constrained connections. PLoS One , 6, e20536. Google Scholar CrossRef Search ADS PubMed  57. Barnett L., Di Paolo E. & Bullock S. ( 2007) Spatially embedded random networks. Phys. Rev. E , 76, 056115. Google Scholar CrossRef Search ADS   58. Kosmidis K., Havlin S. & Bunde A. ( 2008) Structural properties of spatially embedded networks. Europhys. Lett. , 82, 48005. Google Scholar CrossRef Search ADS   59. Robins G., Pattison P., Kalish Y. & Lusher D. ( 2007) An introduction to exponential random graph ($$p*$$) models for social networks. Soc. Netw. , 29, 173– 191. Google Scholar CrossRef Search ADS   60. Wong L. H., Pattison P. & Robins G. ( 2006) A spatial model for social networks. Phys. A Stat. Mech. Appl. , 360, 99– 120. Google Scholar CrossRef Search ADS   61. Haenggi M., Andrews J. G., Baccelli F., Dousse O. & Franceschetti M. ( 2009) Stochastic geometry and random graphs for the analysis and design of wireless networks. IEEE J. Selected Areas Commun. , 27, 1029– 1046. Google Scholar CrossRef Search ADS   62. Brockmann D. & Helbing D. ( 2013) The hidden geometry of complex, network-driven contagion phenomena. Science , 342, 1337– 1342. Google Scholar CrossRef Search ADS PubMed  63. Hoff P. D., Raftery A. E. & Handcock M. S. ( 2002) Latent space approaches to social network analysis. J. Am. Stat. Assoc. , 97, 1090– 1098. Google Scholar CrossRef Search ADS   64. Matias C. & Robin S. ( 2014) Modeling heterogeneity in random graphs through latent space models: a selective review. ESAIM Proc. Surv. , 47, 55– 74. Google Scholar CrossRef Search ADS   65. Serrano M. A., Krioukov D. & Boguná M. ( 2008) Self-similarity of complex networks and hidden metric spaces. Phys. Rev. Lett. , 100, 078701. Google Scholar CrossRef Search ADS PubMed  66. Volz E. M. ( 2008) SIR dynamics in random networks with heterogeneous connectivity. J. Math. Biol. , 56, 293– 310. Google Scholar CrossRef Search ADS PubMed  67. Miller J. C. ( 2014) Epidemics on networks with large initial conditions or changing structure. PLoS One , 9, e101421. Google Scholar CrossRef Search ADS PubMed  68. Fisher R. A. ( 1937) The wave of advance of advantageous genes. Ann. Eugenics , 7, 355– 369. Google Scholar CrossRef Search ADS   69. Kolmogorov A. N., Petrovskii I. G. & Piskunov N. S. ( 1991) A study of the diffusion equation with increase in the amount of substance, and its application to a biological problem. In: Tikhomirov V. M. ed. Selected Works of A. N. Kolmogorov . Dordrecht: Springer. 70. Coville J. & Dupaigne L. ( 2005) Propagation speed of travelling fronts in non local reaction–diffusion equations. Nonlinear Anal. Theory Meth. Appl. , 60, 797– 819. Google Scholar CrossRef Search ADS   71. Berestycki H., Nadin G., Perthame B. & Ryzhik L. ( 2009) The non-local Fisher-KPP equation: travelling waves and steady states. Nonlinearity , 22, 2813. Google Scholar CrossRef Search ADS   72. Hallatschek O. & Fisher D. S. ( 2014) Acceleration of evolutionary spread by long-range dispersal. Proc. Natl. Acad. Sci. U.S.A , 111, E4911– E4919. Google Scholar CrossRef Search ADS PubMed  73. Miller J. C. & Hagberg A. ( 2011) Efficient generation of networks with given expected degrees. In Proceedings of the 8th International Workshop on Algorithms and Models for the Web Graph . pp. 115– 126. 74. Panja D. ( 2004) Effects of fluctuations on propagating fronts. Phys. Rep. , 393, 87– 174. Google Scholar CrossRef Search ADS   75. Penrose M. ( 2003) Random Geometric Graphs , vol. 5. Oxford: Oxford University Press. Google Scholar CrossRef Search ADS   76. Joos F., Perarnau G., Rautenbach D. & Reed B. ( 2016) How to determine if a random graph with a fixed degree sequence has a giant component. IEEE 57th Annual Symposium on Foundations of Computer Science , pp. 695– 703. 77. Bollobás B., Janson S. & Riordan O. ( 2007) Spread-out percolation in d. Random Struct. & Algorithms , 31, 239– 246. Google Scholar CrossRef Search ADS   78. Batagelj V. & Brandes U. ( 2005) Efficient generation of large random networks. Phys. Rev. E , 71, 036113. Google Scholar CrossRef Search ADS   © The authors 2018. Published by Oxford University Press. This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact journals.permissions@oup.com

Journal

Journal of Complex NetworksOxford University Press

Published: Mar 13, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off