Epidemic dynamics on higher-dimensional small world networks.

Dimension governs dynamical processes on networks. The social and technological networks which we encounter in everyday life span a wide range of dimensions, but studies of spreading on finite-dimensional networks are usually restricted to one or two dimensions. To facilitate investigation of the impact of dimension on spreading processes, we define a flexible higher-dimensional small world network model and characterize the dependence of its structural properties on dimension. Subsequently, we derive mean field, pair approximation, intertwined continuous Markov chain and probabilistic discrete Markov chain models of a COVID-19-inspired susceptible-exposed-infected-removed (SEIR) epidemic process with quarantine and isolation strategies, and for each model identify the basic reproduction number R 0 , which determines whether an introduced infinitesimal level of infection in an initially susceptible population will shrink or grow. We apply these four continuous state models, together with discrete state Monte Carlo simulations, to analyse how spreading varies with model parameters. Both network properties and the outcome of Monte Carlo simulations vary substantially with dimension or rewiring rate, but predictions of continuous state models change only slightly. A different trend appears for epidemic model parameters: as these vary, the outcomes of Monte Carlo change less than those of continuous state methods. Furthermore, under a wide range of conditions, the four continuous state approximations present similar deviations from the outcome of Monte Carlo simulations. This bias is usually least when using the pair approximation model, varies only slightly with network size, and decreases with dimension or rewiring rate. Finally, we characterize the discrepancies between Monte Carlo and continuous state models by simultaneously considering network efficiency and network size.

Introduction

Dimension is an important physical property, governing widespread phenomena such as diffusion and modes of vibration [1]. In the context of networks, dimension not only governs dynamical processes [1], [2] and topological properties [3], but even suffices to approximate general structure [4]. While physical objects are constrained to three or fewer dimensions, the dimensions of social, technological, biological and other networks span a wide spectrum [5], [6], [7]. A proper understanding of dynamical network phenomena, such as random walks [8], [9], synchronization and consensus [10], [11], [12], [13], heat conduction [14], [15], information dissemination [16], [17], [18], [19], [20], [21], competition and cooperation [22], [23], [24], [25] and the spread of rumours [26], [27], opinions [28] and diseases [29], [30], [31], [32], [33], [34], [35] including the ongoing COVID-19 pandemic [36], [37], requires assessment of the process across a broad range of dimensions and other network properties. This requires the development and application of flexible higher-dimensional network models.

The Watts-Strogatz small world model [38] captures tension between order and randomness, and is one of the most influential and widely utilized models in network science [39]. The original model was defined by rewiring links of a toroidal lattice of dimension one (i.e., a ring lattice) and arbitrary even mean degree. Rewiring is performed such that, regardless of the rewiring probability, the mean degree of the network remains constant and the degree of each node is at least half the mean degree, making isolated nodes impossible and disconnections unlikely. Since the original small world model was proposed, methods have been developed to generalize to lattices of higher dimension [4] (especially two-dimensional lattices [40], [41], [42], [43]), but each which we have encountered sacrifices at least one important property of the original small world model. Most models applicable to higher dimensions add links instead of rewiring, which changes the mean degree, and/or do not consider the case of mean degree other than twice the network dimension [44], [45], [46], [47]. In this paper we develop a small world model, defined by rewiring links of a toroidal lattice of arbitrary positive integer dimension and arbitrary even degree no smaller than twice the lattice dimension (code is available in Ref. [48]). This model follows the foundational Watts-Strogatz small world model in preserving the mean degree of the original toroidal lattice and ensuring that the degree of each node is at least half the mean degree. We examine the dependence of important structural properties on parameters and implementation choices of the network model.

A range of important phenomena and operations can be represented as dynamical processes on networks [49], including spreading, network synchronization and control of complex systems [50]. However, despite the wide spectrum of dimensions apparent in real networks, dynamical processes in higher-dimensional networks are rarely analyzed. In this contribution we consider susceptible-exposed-infected-removed (SEIR) disease spreading in -dimensional spatial networks, incorporating quarantine and isolation methods. There exist other iconic disease spreading models, such as the susceptible-infected-removed (SIR) or susceptible-infected-susceptible (SIS) models, but we focus on the SEIR model because it is frequently chosen to describe the ongoing COVID-19 pandemic [43], [51]. We use four probabilistic continuous state theories - namely, mean field (MF) [52], pair approximation (PA) [53], intertwined continuous Markov chain (ICMC) [54] and probabilistic discrete Markov chain (PDMC) [55] - to establish analytically tractable spreading dynamical systems. For each continuous state system we deduce the basic reproduction number , which governs whether an infinitesimal level of infection in an initially susceptible population will increase. Moreover, we compare Monte Carlo (MC) and the four continuous state methods, examining how predicted outcomes and discrepancies between methods depend on epidemic, control and network model parameters, and identify systematic bias shared by each of the continuous state methods.

The output of the MF, PA, ICMC and PDMC continuous state models is useful, but just as important is an understanding of when these models are valid and the conditions under which systematic bias occurs [56], [57]. Some discrepancies between MC and continuous state models can be explained by finite size effects [58], [59], in which case errors should decrease with the size of the network. An increased network size allows each MC simulation to involve a larger number of instances of discrete processes, which could indeed bring them closer to continuous state models. However, the ease with which distinct nodes on the network can interact with each other, which is represented by the network efficiency, could also be expected to play a role. To improve understanding of the conditions upon which model accuracy depends, we investigate how network size and network efficiency influence the difference between the final affected ratio predicted by continuous state models and that determined by MC. Furthermore, we formulate a simple network measure, called the network capacity, which is intended to capture effective network size and which we find better describes the discrepancy exhibited by the most accurate continuous state model.

The remainder of the paper is arranged as follows. The flexible higher-dimensional small world model is presented in Section 2, where the structural properties of the resultant small world networks are also examined and where we show that the higher-dimensional small world model is important both for fitting topological properties and predicting dynamical characteristics of social networks. In Section 3 we focus on dynamics, and formulate a susceptible-exposed-infected-recovered (SEIR) epidemic model designed to reflect current understanding of the ongoing COVID-19 pandemic. In particular, this SEIR model includes the possibility for nodes in the exposed (E) state to spread disease. We also derive the MF, PA, ICMC and PDMC continuous state probabilistic descriptions of the process, and identify the basic reproduction number corresponding to each. We use Section 4 to perform numerical simulations illustrating how the spread of disease in higher-dimensional networks depends upon network parameters such as size, mean degree and, especially, dimension. We also investigate how the predictions of different continuous state models compare with MC simulation, highlight a systematic bias which is frequently shared by each continuous state method, and examine how this bias depends on model parameters and network properties. Finally, in Section 5, we conclude.

Higher-dimensional network model and properties

In this section we define the higher-dimensional small world model, investigate the dependence of its properties on model parameters, and show that dimension is important in fitting and predicting the characteristics of real world social networks.

Network measures

First we define relevant network measures. The number of nodes in a network is called size and denoted . The mean network distance between distinct nodes iswhere is the number of edges in the shortest path between nodes and . The clustering coefficient is the probability for a randomly chosen connected ordered triple of nodes to form a triangle:where is the adjacency matrix of the network. Network efficiencyis the mean rate of communication between distinct nodes in the network, assuming that information propagates at a rate of one link per time step. To characterize difference between continuous state and Monte Carlo approaches, we also define the network capacitywhich comprises the mean over all nodes of the total number of communications with other nodes per time step. As the product of the number of nodes () and the mean number of communications per node per time (), the capacity equals the total number of communications in the network per unit time. This formulation is intended to capture the effective network size from the perspective of continuous state approximations (Section 3.2) which rely on the average, over many pairs of nodes, of interactions between nodes.

We also characterize networks by their spectrum; the eigenvalues of the adjacency matrix. The largest eigenvalue of a graph, also called the index or spectral radius, plays an important role in dynamic processes on networks. In disease spreading, the spectral radius directly determines epidemic thresholds [54], [60].

The eigenvalues and also impact dynamic or structural properties. Increasing the difference between the first and second largest eigenvalues, which is called the spectral gap, can improve expansion and connectivity [61], [62].

The smallest eigenvalue provides information about bipartite subgraphs, independence number, chromatic number and the maximum cut [62], [63], [64], [65].

Definition of the higher-dimensional small world model

Next, we define the higher-dimensional small world network model. A small world network of dimension , even mean degree and rewiring probability parameter is obtained as follows. First, we generate a -torus1 of nodes each connected to its nearest neighbours. To do so, we start with nodes distributed in a grid of side length with periodic boundary conditions in each of the directions. Initially, an edge is placed between nodes and if and only if . In this expression, is a version of Euclidean distance2 which respects the periodic boundary conditions in each dimension, in which each dimension is treated equivalently, and is the minimum value for which the resulting network has degree at least . If the degree of this network exceeds the desired degree , then links directed in the clockwise direction3 from a single node are chosen uniformly at random from among the links which span the largest Euclidean distance. These links are removed from each node, which reduces to the degree of the network, while keeping the neighbourhood of each node the same up to translation.

Given a -torus with mean degree , we consider each node in turn, and consider consecutively each of the edges which are directed clockwise starting from node . With probability , we reconnect this edge to a vertex selected from among all nodes with which is not already connected, chosen with probability proportional to , where is the Euclidean distance between the two nodes, and control the balance between uniform random and distance dependent node choice, and the exponent determines how fast probability decays with Euclidean distance. We consider two choices of these parameters: , and arbitrary, which, because in this case rewiring does not depend on distance, we refer to as random rewiring; and , and , which we refer to as distance dependent rewiring.

To choose a small world network of dimension and size about , we generate a small world network of size , where is the integer closest to the real number . In Section 4 we will sample from an ensemble of small world network models which, in particular, requires random choices of network size , rewiring probability and dimension . To randomly choose a small world network with mean degree , size between about and and rewiring rate bounded below by , we proceed as follows. First we randomly choose a rewiring rate such that the rate of choosing a rewiring rate within any subinterval of is proportional to the logarithmic width of that subinterval. Next, we choose a dimension uniformly at random from the set . Subsequently, we choose an approximate network size with probability proportional to ; such that the probability of choosing within an interval of is proportional to the interval’s logarithmic width. The small world network we finally generate has dimension , rewiring probability , and size about (i.e., size ).

Properties of the small world model

Here we examine the properties of small world models of dimension between 1 and 5. Three kinds of networks are considered: (1) lattices; (2) randomly rewired higher-dimensional small world networks; and (3) higher-dimensional small world networks with distance dependent rewiring. In Fig. 1 , we present regular lattice networks, of (mean) degree and dimension from to 5, which are the outcome of the first stage of the algorithm to generate higher-dimensional small world networks. With mean degree fixed at , dimension would lead to a disconnected network with components, each comprising a lattice of dimension having nodes in each direction. For each dimension we choose a network size sufficiently large that the subgraph induced by a node and its neighbours would not change with increasing network size .

Fig. 1

Regular lattices of (mean) degree and dimension (a) , (b) , (c) , (d) and (e) . The subgraph induced by a node and its neighbours is shown with larger, black nodes and thicker, red links. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

We illustrate how the difference between random and distance dependent rewiring depends on network structure via Figs. 2 and 3 . These figures show interpolation between a lattice and a fully rewired network in dimension . It is particularly interesting to compare networks completely rewired () randomly (Fig. 2(c)) and with distance dependence (Fig. 3(c)). The structure resulting from complete distance dependent rewiring is influenced by the original lattice, but this is not evident after complete random rewiring. In Fig. 4 we present the variation of important structural properties of regular lattices of various dimension and (mean) degree . Figure 4(a) showcases global clustering coefficient , network efficiency and mean network distance for a range of dimensions in lattice networks of fixed (mean) degree and network size about . Clustering coefficient decreases steadily as dimension increases. Network efficiency tends to grow as dimension increases, although a slight decrease is evident as dimension increases to - this decrease results from an increase in network size from when to when . Mean network distance decreases rapidly from dimension 1 to 2, but only slowly from dimensions 2 through to 5. In Fig. 4(b) we present the variation of these three network properties with (mean) degree , for fixed dimension and network size . For (mean) degree , the network is disconnected with some infinite internode distances, leading to efficiency . Similarly, for (mean) degree , there is a complete absence of triangles, ensuring that clustering coefficient vanishes. We see that, in the domain of degree where they are non-zero, clustering coefficient increases with growing (mean) degree more slowly than network efficiency .

Fig. 2

Using the small world network with random rewiring to interpolate between a two-dimensional lattice (rewiring rate ) and a random network (). Network size and mean degree .

Fig. 3

Using the small world network with distance dependent rewiring to interpolate between a two-dimensional lattice (rewiring rate ) and a random network (). The probability for a node to be a recipient of a rewired link decreases with the periodic Euclidean lattice distance as (Section 2.2). Network size and mean degree .

Fig. 4

Structural properties of regular lattices networks depend on dimension and mean degree . (a) Variation with dimension of clustering coefficient , network efficiency , and mean network distance , for regular lattices with (mean) degree and about nodes. (b) Variation with (mean) degree of clustering coefficient , network efficiency , and mean network distance , for regular lattices of dimension with nodes.

We reveal the impact of rewiring rate and rewiring strategy (distance dependent or random) in Fig. 5 . Figure 5(a) shows the variation with rewiring rate of the clustering coefficient and network efficiency, in dimension and with mean degree , while Fig. 5(b) shows the variation of mean network distance . Under either rewiring strategy, network efficiency increases with rewiring rate , while clustering coefficient and mean network distance decrease. For a particular rewiring rate , distant-dependent rewiring leads to a higher clustering coefficient and lower efficiency than random rewiring. Similarly, the mean network distance between distinct nodes is lower under random rewiring than distance-dependent rewiring. Results for other dimensions are similar. In Fig. 5(c,d) we show how the largest, second largest and smallest eigenvalues, , and , vary with random and distance dependent rewiring rate in dimension to 5 with mean degree . The largest eigenvalue of a graph is at least as large as the average degree [66] and, for regular graphs, exactly equals the (mean) degree [62]. Accordingly, the largest eigenvalue is in the absence of rewiring (), but grows as rewiring rate increases.

Fig. 5

Network structural properties depend on rewiring rate and rewiring strategy. For a network with dimension , mean degree and nodes, variation with random and distance dependent rewiring rate of: (a) clustering coefficient and network efficiency ; and (b) mean network distance . For networks of five different dimensions, size about and mean degree , the largest eigenvalue , second largest eigenvalue and negative of the smallest eigenvalue under: (c) random rewiring; and (d) distance dependent rewiring. The values of , and shown comprise the mean over twenty independently generated graphs.

For fixed rewiring rate , as dimension increases, the spectral gap increases, indicating better connectivity in higher dimension. The magnitude of the smallest eigenvalue also increases with rewiring rate and, under random rewiring, and converge.

As rewiring rate grows, the second largest and smallest eigenvalues, and , decrease more slowly under distance dependent rewiring than random rewiring. Under distance dependent rewiring, the second smallest eigenvalue even increases in dimension , and initially also increases in distance , maintaining a small spectral gap and, hence, high diameter and low connectivity. In particular, the dependence on rewiring rate (of clustering coefficient , efficiency , mean distance and the eigenvalues , and ) shows that, for the same fraction of rewired links, distance dependent rewiring can provide small world networks with structural properties closer to those of regular lattices.

Utility of the higher-dimensional small world model

In this section we show that the higher-dimensional small world model with random rewiring improves fits and predictions of properties of observed real world networks. We focus on the measures considered in Ref. [38], which proposed the original one-dimensional small world network model.4 Given an observed network, first, we fit the parameters of a small world model to the topological properties of clustering coefficient and mean network distance . Next, we use the fitted small world model parameters to calculate the dynamical properties examined in Ref. [38], which are defined in terms of a deliberately simplified susceptible-infected-removed (SIR) model. In this SIR model, spreading begins from a single, randomly chosen infected node in an otherwise completely susceptible network. In each time step, each infected node infects each of its susceptible neighbours with probability and, simultaneously, all infected nodes enter the removed state. The dynamical properties investigated in Ref. [38] are the time required for a maximally infectious disease, with , to spread throughout the entire network, and the minimum rate of infection for which the disease infects at least half the network. Each of these properties is calculated for the observed real world network by taking the mean over 100 independent trials.

The small world parameters we must fit are mean degree , dimension , network size , and rewiring rate . Mean degree is chosen as the even integer , where is the mean degree of the observed real world network and represents the integer closest to . Dimension can then be chosen in the range . Given a choice of dimension , network size is fitted as the power of closest to the observed network size . The rewiring rate is then chosen to minimize the root mean squared relative error in the fit:whereis the relative error in a statistic and is the value of the statistic observed in the real world network. This minimization is performed using the golden section method with variable in the interval , which involves the simplifying assumption that has a single local minimum on the domain considered. For each value of , the objective function is estimated as the mean over 100 independent realisations of the small world network model. When we consider the original, one-dimensional small world model, dimension is fixed at . In contrast, when we consider the higher-dimensional small world model, each dimension is considered, with the dimension finally chosen that which provides the smallest value of .

At this point we have fitted the parameters of the small world model to the observed real world network. We now recalculate topological and dynamical properties , , and as the mean over 100 trials, each from an independently generated small world network with the identified model parameters. The root mean squared relative error to the fitted (topological) properties is recalculated from Eq. (1), while the root mean squared relative error to the predicted (dynamical) properties is calculated as:

We apply this approach to the largest connected components of eighty-four high school friendship networks [68]. Table 1 presents basic properties of these networks, as well as model parameters estimated using the original small world model (fixed dimension ) and the higher-dimensional extension (fitted dimension ). The fitted network dimension is consistently three or less, and has median . The table shows that allowing dimension to exceed one allows network properties to be fit with lower rewiring rate .

Table 1

Properties of the largest connected components of eighty-four US high school friendship networks [68]. From top bottom: observed network size, ; observed mean degree, ; observed clustering coefficient, ; observed mean network distance, ; fitted dimension, ; fitted rewiring rate using dimension (i.e., using original small world model), ; fitted rewiring rate using fitted dimension, (i.e., using the higher-dimensional small world model), .

Property	Minimum	First quartile	Median	Third quartile	Maximum
N^o	25	222.5	387	621	1974
〈k〉^o	3.7	6.2	7.4	8.1	12.0
ϕ^o	0.13	0.17	0.21	0.26	0.67
〈d〉^o	1.5	3.3	3.8	4.2	6.6
D^	1	1	2	2	3
p|_D=1	0.012	0.24	0.29	0.34	0.42
p|_D=D^	0.0014	0.12	0.19	0.25	0.34

In Fig. 6 we compare fitting error (to topological properties) and prediction error (of dynamical properties) using the established, one-dimensional small world model and the developed higher-dimensional small world model. The legends summarize results by stating the number of times each approach is optimal (Opt.), providing the (possibly equal) lowest error. As Fig. 6(a) shows, fitting dimension as well as rewiring rate reduces fitting error in most cases, and rarely increases error. There is a single counter-intuitive case in which fitting error for exceeds that for . This corresponds to a situation in which statistical fluctuations cause dimension to minimize using the mean over 100 trials calculated during the fitting procedure, but this fitted dimension did not minimize during its final evaluation using the mean over 100 new trials. In the majority of cases, as Fig. 6(b) reveals, fitting dimension also reduces error in the prediction of dynamical properties. The networks for which the original small-dimensional model provides (possibly equal) lowest and tend to coincide: for twenty-two networks provides both lowest fitting error and lowest prediction error .

Fig. 6

Access to higher dimensions improves fitting to topological properties and prediction of dynamical properties. Errors in fitting and prediction of properties of eighty-four US high school friendship networks (Table 1) [68]. (a) Error in fitted, topological properties, calculated using Eq. (1). (b) Error in predicted, dynamical properties calculated using Eq. (2). The original small world network model is restricted to dimension , but the higher-dimensional network model allows the use of a the fitted dimension . The number of times each approach is optimal (Opt.) and yields the (possibly equal) lowest error is shown in the legends. For readability, the observed real worlds networks are ordered, separately for (a) and (b), by increasing error in the case.

SEIR epidemic model and theoretical analysis

In this section, we establish a COVID-19-inspired SEIR disease spreading model incorporating quarantine and isolation control strategies. We apply four continuous state methods - MF, PA, ICMC and PDMC methods - to analyze SEIR processes in higher-dimensional networks generated using the developed small world model. For each method we determine the basic reproduction number .

SEIR epidemic model

In this section we outline a SEIR model designed to reflect known properties of the ongoing COVID-19 pandemic, as well as the intervention strategies applied to its mitigation. Specifically, the strategies we consider are two of the most effective methods for containing COVID-19: isolation of confirmed infections (based on testing); and quarantine of those who have been exposed to infection (based on contact tracing) [69], [70].

A COVID-19 carrier can propagate the disease to others without experiencing symptoms themselves [71], [72], and person-to-person transmission is possible even during the incubation period [73]. Accordingly, in contrast to many SEIR models in the literature, we allow exposed (E) nodes to have a non-zero probability of spreading disease.

We establish the SEIR disease spreading model incorporating quarantine and isolation strategies, as shown in Fig. 7 . Under the SEIR model, each individual is in one of four states: individuals susceptible (S) to the disease; individuals who have been exposed (E) to the disease and who have not yet developed symptoms but are in an incubation period; infected (I) individuals; and individuals who have recovered from the disease and developed immunity, or for other reasons have been removed (R) from further consideration in the spreading system. In contrast to many established SEIR epidemic models, under our model, nodes are only in an E or I state when their movements are not controlled through quarantine or isolation measures which would prevent them from infecting others. The model is specified by rate parameters which govern the probability of transition between states. The transition from S to E can occur in two distinct ways: S contacting E neighbours and changing into E state with rate ; and S contacting I neighbours and changing into E state with rate . The disease evolves such that individuals in the E state enter the I state with rate (the disease progression rate). Simultaneously, the quarantine strategy allow individuals in state E to recover or be removed, and enter state R without involvement in new infections, with rate . The transition from I to R can also occur in two different ways. Infected nodes are isolated with rate , and in this case enter state R without the opportunity to infect new nodes. Infected nodes can also enter state R by recovery, with rate .

Fig. 7

Structure of the SEIR disease spreading model. : infection rate from E nodes; : infection rate from I nodes; : disease progression rate; : quarantine rate; : recovery rate; : isolation rate.

As much as possible, we use rate parameters relevant to the COVID-19 case, although these are not known precisely, and vary with time and location [71], [74], [75], [76]. A seven day incubation period of COVID-19 [43], [77], [78] leads to a progression rate , and a fourteen day convalescence period [43], [79] provides a recovery rate . For other parameters we consider a range of reasonable values.

Theoretical analysis

In this section we apply four distinct continuous state methods to describe the SEIR epidemic spreading process and identify the basic reproduction number .

Mean field

Higher-dimensional small world networks have Poisson-like degree distributions. Therefore, we begin by considering a continuous state method developed for homogeneous uncorrelated networks. The MF formulation of the dynamical spreading system iswhere , , and denote the fraction of individuals in the S, E, I and R state respectively, so that .

We will use the next generation matrix (NGM) [80] method to deduce the basic reproduction number. Let and be the subpopulations in each of the possessive compartment (the states, E and I, which possess the disease) and non-possessive compartment (the states, S and R, which do not possess the disease) respectively:

The compartmental model is represented by the vector functionswhere represents the rate of transmissions (from non-possessive to possessive compartments) and describes transitions (between possessive or from possessive to non-possessive), such that . The SEIR model has a disease free equilibrium (DFE) at which the population remains in the absence of disease, and this equilibrium solution has the form . When the infectious disease is first introduced, at most a small number of individuals will be in the removed state, and so we can use the approximation at the DFE. It follows that the Jacobian matrices of and at the DFE are

The basic reproduction number of model (3) can be expressed in terms of spectral radius aswhere denotes the spectral radius of the square matrix .

Pair approximation

The standard MF Eq. (3) assume that the states in a neighborhood are independent, i.e. that there are no dynamical correlations. The MF method considers each node in isolation (i.e., first order interactions) but more accurate descriptions are possible by incorporating higher order interactions, where the order of an interaction is the number of nodes which it involves. A technique incorporating higher order interactions is the moment closure approximation [81], [82]. When quantities summarising third-order interactions are approximated in terms of quantities describing interaction of order two or lower (system equations are “closed” at the level of pairs) the moment closure method can be called PA [82]. The PA and MF methods both assume randomly mixing populations, with the distinction that PA assumes that pairs of individuals, rather than the individuals themselves, randomly mix [83]. The PA method is most appropriate for networks with little variance in the degrees, i.e., for which the number of neighbours of individuals do not differ too much from the average number of neighbours per individual [84].

Here, we extend the PA method to our SEIR model, adopting the following notation: let be the number of sites that are of type X, let be the number of ordered X-Y pairs and, similarly, let be the number of ordered X-Y-Z triples. The following differential equations describe the evolution of first order (node level) quantities according to the PA description of the SEIR model:The corresponding differential equations for second order quantities (describing states of pairs of nodes) areIn order to close Eqs. (4) and (5), we must express the third-order quantities, like , etc., in terms of the second-order state variables . We employ the homogeneous network approximation [81], [83] The PA system comprises the thirteen equations given by systems (4) and (5), together with the preceding approximation of third order interactions.

The basic reproduction number can be estimated by introducing correlations between individuals of distinct type and :This definition allows system (4) to be rewrittenBy applying the NGM method as in Section 3.2.1, we can estimate the basic reproduction number asAssuming that exposed and infected individuals are chosen randomly at the initial time, . Therefore, under PA we estimate the same basic reproduction number as under the MF method.

Intertwined continuous Markov chain

In this section we present the ICMC [54] method to analyze disease spreading in networks. Let and represent the probability at time for the node to be S, E, I and R respectively. These probabilities satisfy the normalisation condition .

The probabilities and for node in state S not to be affected by any E and I neighbours are respectivelyThe ICMC system can therefore be expressed

We use the NGM method to deduce the basic reproduction number . Analogously to Section 3.2.1, and , are the possessive compartments, and the compartmental model can be writtenwhere and . Analogously to Section 3.2.1, the disease free equilibrium (DFE) solution of ICMC method is , with and . The Jacobian matrices of and at DFE point are thenwhere is the identity matrix. The basic reproduction number of system (5) is , or

Since spectral radius equals mean degree in the absence of rewiring () and increases only slightly with rewiring rate (Fig. 5(c,d)), the value of the basic reproduction number estimated via the ICMC method is similar to the result obtained using MF.

Probabilistic discrete Markov chain

The fourth and final continuous state approach which we consider is the PDMC [55] method. Let denote the probability that node is in state X at time , and observe that . The PDMC update rules for node arewhere are given by Eq. (7).

To identify under the PDMC method, following Ref. [55], [85], we consider introducing into an otherwise purely susceptible system small levels of exposed and/or infected nodes. We then identify the critical epidemiological parameters and network properties for which an approximately stable solution exists. Conditions more (less) favourable for epidemics than this critical point will correspond to levels of disease which initially grow (shrink). Accordingly, near the critical point, the probability for any node to be in states E and I are , and . By Eq. (7),Inserting this into system (9), we obtainBy eliminating we reachwhereThus the solution of Eq. (10) reduces to an eigenvalue problem for the matrix , whose elements are . A nontrivial solution to Eq. (10) requires the largest eigenvalue of to equal one, which allows us to identify the basic reproduction number as the spectral radius of ; . This result for is the same as under the ICMC method.

Numerical simulation

In this section, we present numerical simulations illustrating the SEIR epidemic spreading model in -dimension small world networks. Unless stated otherwise, we consider networks of size about and mean degree , and use epidemic parameters , , , and , . Initial conditions are chosen as . Under ICMC and PDMC, for greater realism, we randomly seed the predominately susceptible system with the specified number of exposed, infected and recovered individuals such that each node initially is entirely in a single state; , . For PA, we assume initial conditions corresponding to perfect mixing: , where . Unless indicated otherwise, results for MC comprise the mean over fifty independent trials.

In Fig. 8 we show how the disease evolves in time according to our five methods of description, MF, PA, ICMC, PDMC and MC, in a network of size , dimension , mean degree and random rewiring rate . Figure 8(a,b) show the change with time of the infected and exposed node fraction and . These fractions are lower under MC than under other methods, and PA has the least extreme deviations from MC. Figure 8(c,d) show and ; is the fraction of pairs of (ordered) pairs of connected nodes which comprise a node in state X and a node in the state Y, where, to calculate link fractions under MF, we assume perfect mixing. In Fig. 8(e,f) we show the variation with time of the correlation , for each of the five methods. The fractions of pairs of SE and SI, i.e., and , and the SE and SI correlations and are lowest for the MC method, and PA has the least extreme deviations from MC. Together, these results show that the continuous state methods MF, PA, ICMC and PDMC systematically overestimate levels of disease compared with MC, but that PA offers the least extreme differences.

Fig. 8

Different descriptions of the SEIR spreading process predict distinct time evolution. Variation with time according to five different methods: MF, PA, ICMC, PDMC and MC. (a) Exposed node fraction ; (b) infected node fraction ; (c) SE pair fraction ; (d) SI pair fraction ; (e) SE correlation ; and (f) SI correlation . Spreading takes place on a small world network with size , dimension , mean degree and random rewiring rate .

Dependence on epidemic model parameters

We use Fig. 9 to show the final number of people affected by the epidemic, which is represented by the final affected ratio , depends on infection rates and for a network of size , mean degree and dimension , with random rewiring rate .

Fig. 9

The final number affected by the disease increases with infection rate. (a)–(e) Variation with infection rates and of the final affection ratio when using the (a) MF; (b) PA; (c) ICMC; (d) PDMC; and (e) MC methods. Level sets of the final affected density are shown with black dashed lines, while the white dotted line highlights parameters for which the basic reproduction number is unity. (f) Variation with of the final affection ratio , with . Spreading takes place on a small world network with dimension , size , mean degree and random rewiring rate . Results for Monte Carlo comprise the mean over twenty independent trials.

Under any method, the final affected ratio increases as either infection rate increases. Comparing Fig. 9(a)–(e) shows that, across the spectrum of infection rates considered, the MF, ICMC and PDMC methods lead to similar final affected ratios , which are higher than those under MC, while PA leads to outcomes intermediate between MC and the other continuous state methods. This pattern is highlighted by Fig. 9(f), which shows how the final affected ratio varies with infection rate for each of the five methods considered, for fixed infection rate . Figure 9(a)–(d) also depict the critical values of infection rates around which the basic reproduction number under each continuous state method achieves its critical value . This threshold is similar under any of the continuous state methods and, for the MF, ICMC and PDMC methods, is roughly parallel to the level sets of the final affected ratio .

In Fig. 10 we show how the final affected ratio depends upon epidemic intervention parameters. Figure 10(a) shows how depends on quarantine rate , while Fig. 10(b) shows the variation with isolation rate . Each of the five modelling methods predicts rapid decreases in the final affected ratio as either the quarantine rate or isolation rate increases. The figure shows that quarantine and isolation strategies are effective in decreasing the final number of affected individuals. The figure also shows that the four methods MF, PA, ICMC and PDMC overestimate the final affected ratio compared with MC, but that PA provides a more accurate result.

Fig. 10

The final number affected by the disease decreases with quarantine rate and isolation rate . Variation of final affected ratio with (a) quarantine rate and (b) isolation rate for a small world network of dimension , size , mean degree and random rewiring rate .

The four continuous state methods predict far stronger dependencies on epidemic parameters than are observed under MC. This is especially evident for infection rates and , but is also true of control parameters and . The relative insensitivity of MC to changes in epidemic parameters may arise because continuous state approximations allow arbitrarily small fractions of infections to continue propagating but, under MC, any level of infection which is not one is zero. This would tend to attenuate increases in final affected ratio observed under MC as epidemic parameters are changed. Each of the continuous state models we consider has so far tended to overestimate the final affected ratio , although the PA method does so least.

Dependence on network model parameters

Now we examine the dependence of the final affected ratio on parameters of the small world network model. Figure 11 presents the variation of final affected ratio with mean degree for five dimensions from to 5. The figure reveals that mean degree strongly affects final outcomes. Specifically, Fig. 11(a)–(e) reveal that, in any dimension, the final affected ratio grows with mean degree . Moreover, we notice that MF, ICMC and PDMC provide quite similar results, and under widespread conditions predict higher final affected ratios than MC. The PA method usually overestimates final affected ratio compared with MC, but provides results closer to MC than are provided by the other three continuous state methods. In particular, in dimension or 5, the final outcomes are similar under PA and MC. Each of the continuous state methods is insensitive to dimension, but under MC, the final affected ratio grows with increasing dimension . Thus, analyses based on MC simulations with the original, one-dimensional, small world model have the potential to underestimate the severity of an epidemic.

Fig. 11

The final number affected by the disease increases with mean degree. (a-e) Variation with mean degree of the final affected ratio according to the MF, PA, ICMC, PDMC and MC methods for networks of size about , randomly rewired with probability , and with dimension (a) , (b) , (c) , (d) , and (e) . (f) The difference between MC and other methods for the network of dimension .

In Fig. 11(f) we highlight a trend present throughout Fig. 11(a–e). Figure 11(f) shows that the difference between MC and the continuous state methods is not monotonic in mean degree . Instead, the discrepancy initially increases with mean degree , reaches a peak around , and decreases again. The smaller discrepancy for low degree may arise because interactions are more local and easier to approximate. In turn, the smaller deviation for higher degree can be explained because the network is closer to complete, which corresponds better to the assumptions underlying mean field theory. The figure also illustrates that, in dimension , as for other dimensions, the difference from MC is less for PA than the other three continuous state methods, MF, ICMC and PDMC.

In Fig. 12 we illustrate the effect of random or distance dependent rewiring on epidemic spreading in higher-dimensional networks, focusing on mean degree and networks of size about . For all five methods in dimension , Fig. 12(a) and (b) illustrate the variation of the final affected ratio with random and distance dependent rewiring rate respectively. Under MC we observe that the final affected ratio grows with increasing rewiring rate but, when using the continuous state methods MF, PA, ICMC and PDMC, outcomes vary little with rewiring rate . For large random or distance rewiring rate , the final affected ratio under MC approaches that under PA but, for the other three continuous state methods, a difference persists.

Fig. 12

The final number affected by the disease increases with rewiring rate , but this variation decreases with dimension . Change with rewiring rate of affected node fraction for networks with mean degree and nodes: (a,b) in dimension using MF, PA, ICMC, PDMC and MC methods; and (c,d) using MC in dimensions to 5. Rewiring is: (a,c) random; and (b,d) distance dependent.

Figure 12 (c,d) shows how the final affected ratio changes with random and distance dependent rewiring rate under MC and in dimensions to 5 (under the four continuous state methods, varied little with rewiring rate in any dimension). For fixed rewiring rate , the final affected ratio grows as dimension increases from , the dimension of the original small world model. For lower dimensions , the final affected ratio grows as rewiring rate increases. When the random rewiring rate is high, the final number of affected nodes varies little with dimension , but under distance dependent rewiring even large rewiring rates allow substantial disparities between dimensions. Networks of similar size can accommodate larger distances in smaller dimension, which makes more significant the decision to rewire in a distance-dependent way. For lower dimension, whether the epidemic affects the majority of nodes of a fully rewired network (for the parameters which we consider) hinges on whether rewiring is performed randomly or in a distance dependent manner. For higher dimension, this difference is less obvious. Thus, Fig. 12(c,d) illustrate that both dimension and the nature of shortcuts impact epidemic outcomes, and that long-distance connections increase spreading.

In Fig. 13 we investigate the dependence on network size of the final affected ratio for a diverse selection of dimension (), mean degree () and random rewiring rate (). Here, we use the initial conditions , . For each of the five methods considered, results vary little as network size changes. In particular, MF, ICMC and PDMC consistently yield similar results which comprise an overestimate of the final affected ratio according to MC, while PA is usually intermediate between MC and the other continuous state methods (Fig. 13(b,c)) or, under some conditions, slightly underestimates (Fig. 13(a)). In particular, the figure suggests that the discrepancies between methods are not due to finite size effects.

Fig. 13

The final number affected by the disease varies little with network size . Variation with network size of affected node fraction using MF, PA, ICMC, PDMC and MC methods, for networks with dimension, mean degree and random rewiring rate: (a) , and ; (b) , and ; and (c) , and . .

In Section 2 we showed that important network properties vary substantially with rewiring rate and dimension , but this section shows that these substantial changes in topological properties have limited effect on the outcome of continuous state models. This insensitivity is inconsistent with MC simulations, the results of which vary considerably with network model parameters. We might observe this difference because under MC, the disease cannot propagate via fractional levels of infection - any level of infection less than unity must be zero. This makes MC simulations more dependent on the ability of individual infected or exposed nodes to continue spreading, and this ability will likely hinge on topological properties. As previously observed when varying epidemic model parameters (Section 4.1), each of the continuous state models tends to overestimate the final affected ratio , but the PA method does so least.

Characterization of bias

We have illustrated how, under a wide range of conditions, the four continuous state approximations which we consider - MF, PA, ICMC and PDMC - exhibit similar biases relative to the discrete MC approach. We will now investigate how such discrepancies in final affected density can be characterized using network properties. To do so, we consider randomly chosen small world networks of fixed mean degree , size between about and and rewiring rate bounded below by (Section 2.2), and compare the final affected ratio for each continuous state method X with the mean value of final affected ratio over five MC trials. We consider 100 independently generated networks per mean degree , and use the Kendall correlation coefficient to assess the statistical significance of variations of the magnitude of the deviation with changing network size , efficiency and capacity (Section 2.1).

In Fig. 14 we demonstrate that, for fixed mean degree , the capacity governs the deviation from MC of the most accurate method, PA, more consistently than either the network size or efficiency . As Fig. 14(a) shows, for fixed mean degree , the decrease of error with capacity is significant at the level. Fig. 14(b) does reveal that, as mean degree varies from 2 to 16, it is only rarely the case that the Kendall correlation corresponding to the capacity has higher magnitude than the correlations corresponding to both network size and efficiency , but only capacity consistently leads to correlation which is significant at the level (Fig. 14(c)). We have just noted that, for fixed mean degree , the deviation of PA from MC decreases with capacity . However, Fig. 14(d) shows that, conversely, when different mean degrees are grouped, the deviation instead increases with capacity with Kendall correlation significant at the level. That is, the relationship between capacity and the deviation between PA and MC, grouped by mean degree, provides an example of Simpson’s paradox, in which a trend which appears in each group of data disappears or reverses when the groups are combined.

Fig. 14

For fixed (combined) mean degree , the discrepancy between PA and MC decreases (increases) with increasing capacity . Variation with different network properties of the absolute value of the difference between the mean affected ratio for PA and MC methods for random small world networks of size between about and and rewiring rate bounded below by . (a) For fixed mean degree , variation of the absolute value of the difference versus capacity . (b) Kendall correlation between absolute deviation and network size , efficiency and capacity respectively versus mean degree and (c) levels of significance for hypothesis that error decreases with these network properties. (d) Absolute value of the difference between the mean affected ratio for PA and MC methods versus network capacity , with all mean degrees combined.

Figure 15 pertains to differences in the final affected ratios according to MF and MC (the deviations of ICMC and PDMC from MC are similar). In Fig. 15 we show that, for fixed mean degree , both the efficiency and capacity govern deviation from MC better than the network size for fixed mean degree . We show in Fig. 15(a) that, for fixed mean degree , the decrease of error with efficiency is significant at the level. As shown in Fig. 15(a), for mean degree varying from to 16, in all but one case the Kendall correlation corresponding to the efficiency has the highest magnitude, followed by the correlation corresponding to the capacity , and then the correlation corresponding to network size . However, the Kendall correlation is significant at the level for the efficiency for precisely the same set of mean degrees for which it is significant for the capacity (Fig. 15(c)). The relationship between efficiency and deviation between MF and MC, grouped by mean degree, provides another example of Simpson’s paradox: we see from Fig. 15(d) that when different mean degrees are grouped, the increase with capacity of the deviation is significant at the level.

Fig. 15

For fixed (combined) mean degree , the discrepancy between MF and MC decreases (increases) with increasing efficiency . Variation with different network properties of the absolute value of the difference between the mean affected ratio for MF and MC methods for random small world networks with size between about and and rewiring rate bounded below by (Section 2.2). (a) For fixed mean degree , variation of the absolute value of the difference versus efficiency . (b) Kendall correlation between absolute deviation and network size , efficiency and capacity respectively versus mean degree and (c) levels of significance for the hypothesis that error decreases with these network properties. (d) Absolute value of the difference between the mean affected ratio for MF and MC methods versus network efficiency , with all mean degrees combined.

Conclusion

In this paper we formulated a flexible higher-dimensional small world network model which allowed us to characterize how the spread of disease, represented using the SEIR epidemic model, depends upon network properties including dimension and mean degree. The developed higher-dimensional small world model is defined by rewiring links of a toroidal lattice of arbitrary positive integer dimension and arbitrary even degree no smaller than twice the lattice dimension (code is available in Ref. [48]), and emulates the important properties of the foundational Watts-Strogatz small world model. We investigate structural properties of higher-dimensional graphs by considering three kinds of networks: (1) original lattice networks; (2) higher-dimensional small world networks with random rewiring (independent of distance); and (3) higher-dimensional small world networks with rewiring probability decreasing with increasing distance. We found that important network properties depend on dimension, especially as dimension increases from one to two, with higher dimension leading to lower mean network distance and clustering coefficient but higher network efficiency. Moreover, distance dependent rewiring led to higher clustering coefficient and mean network distance. As random rewiring rate or dimension increased, the largest eigenvalue increased slightly, the second largest eigenvalue decreased, the magnitude of the smallest eigenvalue increased and the spectral gap decreased. Under distance-dependent rewiring, these spectral properties changed more slowly with rewiring rate. In addition, we showed that incorporating dimension into the small world network model improves fits to the topological properties of real networks, as well as improving prediction of their dynamical properties.

To study epidemic spreading dynamics in higher-dimensional small world networks, and incorporating known properties of the COVID-19 epidemic, we extended the SEIR epidemic spreading model while considering quarantine and isolation disease control strategies. Specifically, we applied the continuous state MF, PA, ICMC and PDMC methods, and for each obtained the basic reproduction number . The dependence of on quarantine and isolation strategies parameters indicates the effectiveness of these strategies for controlling outbreaks. We verified the effectiveness of these control strategies with numerical simulations involving the four continuous state methods as well as MC sampling, and in this way corroborate other recent studies [70], [74], [86]. Additional results included the following. First, according to MC simulations, the number affected by the disease increases with dimension. This suggests that analyses based on the original small world model, which was one-dimensional, could underestimate the true impact of an epidemic. Second, for fixed mean degree, the difference from MC is similar for each of the four continuous state methods considered, but usually least for PA, and decreases as dimension increases.

Third, the final affected ratio varies little with network size, either for MC or the four continuous state methods. Fourth, the final affected ratio changes more under continuous state methods than MC simulations as epidemic parameters vary. In contrast, as the dimension or rewiring rate change, the outcomes of MC simulations change more than the predictions of continuous state methods. Fifth and finally, for fixed mean degree, the network efficiency governs the deviation of MF from MC more consistently than the network size , but a new measure - the network capacity - governs most consistently the deviation of PA from MC.

The preceding investigation of disease spreading on higher-dimensional small world networks affects both future epidemiological practice and emerging network science research. We have illustrated the impact of epidemic and control parameters across networks of realistically broad degree and dimension, and this can inform disease control and prevention policies. Just as important for decision-makers is information about the limitations of popular continuous state models, and how these limitations depend on topological characteristics. Finally, the developed flexible higher-dimensional small world model opens the door to a plethora of studies gauging the impact of dimension, mean degree and randomness on structural properties and dynamical processes.

Code and data availability

Code to produce higher-dimensional small world networks and the data which support the findings of this study are openly available at https://github.com/JackMurdochMoore/small-world.