key: cord-0832579-4zaryci1
authors: Wang, Haiying; Moore, Jack Murdoch; Small, Michael; Wang, Jun; Yang, Huijie; Gu, Changgui
title: Epidemic dynamics on higher-dimensional small world networks
date: 2022-05-15
journal: Appl Math Comput
DOI: 10.1016/j.amc.2021.126911
sha: 6aa5204a72e6c17d7f0778d110a849c94f86fff0
doc_id: 832579
cord_uid: 4zaryci1

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 [Formula: see text] , 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.

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 D -dimensional spatial networks, incorporating quarantine and isolation methods. There exist other iconic disease spreading models, such as the susceptible-infectedremoved (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 R 0 , 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 R 0 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.

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.

First we define relevant network measures. The number of nodes in a network is called size and denoted N. The mean network distance d between distinct nodes is

where d i j is the number of edges in the shortest path between nodes i and j. The clustering coefficient φ is the probability for a randomly chosen connected ordered triple of nodes to form a triangle:

where A = a i j is the adjacency matrix of the network. Network efficiency

is 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 capacity χ = Nη, which 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 ( N) 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 1 ≥ 2 ≥ . . . ≥ N of the adjacency matrix. The largest eigenvalue 1 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 1 directly determines epidemic thresholds [54, 60] .

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

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

Next, we define the higher-dimensional small world network model. A small world network of dimension D , even mean degree k and rewiring probability parameter p ∈ [0 , 1] is obtained as follows. First, we generate a D -torus 1 of N = L D nodes each connected to its k nearest neighbours. To do so, we start with L D nodes distributed in a grid of side length L with periodic boundary conditions in each of the D directions. Initially, an edge is placed between nodes i and j if and only if 0 < d 2 (i, j) < r. In this expression, d 2 is a version of Euclidean distance 2 which respects the periodic boundary conditions in each dimension, in which each dimension is treated equivalently, and r > 0 is the minimum value for which the resulting network has degree at least k . If the degree k 0 of this network exceeds the desired degree k , then 1 2 ( k 0 − k ) links directed in the clockwise direction 3 from a single node are chosen uniformly at random from among the links which span the largest Euclidean distance. These 1 2 ( k 0 − k ) links are removed from each node, which reduces to k the degree of the network, while keeping the neighbourhood of each node the same up to translation.

Given a D -torus with mean degree k , we consider each node i in turn, and consider consecutively each of the k/ 2 edges which are directed clockwise starting from node i . With probability p, we reconnect this edge to a vertex selected from among all nodes with which i is not already connected, chosen with probability proportional to κ + ω d 2

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: κ = 1 , ω = 0 and σ arbitrary, which, because in this case rewiring does not depend on distance, we refer to as random rewiring; and κ = 0 , ω = 1 and σ = 3 , which we refer to as distance dependent rewiring.

To choose a small world network of dimension D and size about ˜ N , we generate a small world network of size ˜

is the integer closest to the real number x . In Section 4 we will sample from an ensemble of small world network models which, in particular, requires random choices of network size N, rewiring probability p and dimension D . To randomly choose a small world network with mean degree k , size between about ˜ N min and ˜ N max and rewiring rate bounded below by p min ∈ (0 , 1] , we proceed as follows. 

Here we examine the properties of small world models of dimension D 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 k = 10 and dimension from D = 1 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 k = 10 , dimension D > 5 would lead to a disconnected network with L D −5 components, each comprising a lattice of dimension D = 5 having L = N 1 /D nodes in each direction. For each dimension we choose a network size N sufficiently large that the subgraph induced by a node and its neighbours would not change with increasing network size N.

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 D = 2 . It is particularly interesting to compare networks completely rewired ( p = 1 ) 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 D and (mean) degree k . Figure 4 (a) showcases global clustering coefficient φ, network efficiency η and mean network distance d for a range of dimensions D in lattice networks of fixed (mean) degree k = 10 and network size about N = 3 , 0 0 0 . Clustering coefficient φ decreases steadily as dimension D increases. Network efficiency η tends to grow as dimension increases, although a slight decrease is evident as dimension increases to D = 5 -this decrease results from an increase in network size from N = 2401 = 7 4 when D = 4 to N = 3025 = 5 5 when D = 5 . Mean network distance d 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 k , for fixed dimension D = 2 and network size N = 3 , 025 . For (mean) degree k < 2 D , the network is disconnected with some infinite internode distances, leading to efficiency η = 0 . Similarly, for (mean) degree k ≤ 2 D , 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 k more slowly than network efficiency η.

We reveal the impact of rewiring rate p and rewiring strategy (distance dependent or random) in Fig. 5 . Figure 5 (a) shows the variation with rewiring rate p of the clustering coefficient and network efficiency, in dimension D = 2 and with mean degree k = 10 , while Fig. 5 (b) shows the variation of mean network distance d . Under either rewiring strategy, network efficiency η increases with rewiring rate p, while clustering coefficient φ and mean network distance d decrease. For a particular rewiring rate p, distant-dependent rewiring leads to a higher clustering coefficient φ and lower efficiency η than random rewiring. Similarly, the mean network distance d 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, 1 , 2 and N , vary with random and distance dependent rewiring rate p in dimension 

Using the small world network with random rewiring to interpolate between a two-dimensional lattice (rewiring rate p = 0 ) and a random network ( p = 1 ). Network size N = 100 and mean degree k = 4 . Fig. 3 . Using the small world network with distance dependent rewiring to interpolate between a two-dimensional lattice (rewiring rate p = 0 ) and a random network ( p = 1 ). The probability for a node to be a recipient of a rewired link decreases with the periodic Euclidean lattice distance d 2 as d 2 D = 1 to 5 with mean degree k = 10 . 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 1 = k = 10 in the absence of rewiring ( p = 0 ), but grows as rewiring rate p increases.

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

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

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 d . 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 λ ∈ [0 , 1] and, simultaneously, all infected nodes enter the removed state. The dynamical properties investigated in Ref. [38] are the time T required for a maximally infectious disease, with λ = 1 , to spread throughout the entire network, and the minimum rate of infection λ Half 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 k , dimension D , network size N, and rewiring rate p. Mean degree k is chosen as the even integer 2 k o / 2 , where k o is the mean degree of the observed real world network and 

is the relative error in a statistic s and s o is the value of the statistic observed in the real world network. This minimization is performed using the golden section method with variable log p in the interval [ − log k N, 0 ] , which involves the simplifying assumption that F has a single local minimum on the domain considered. For each value of p, the objective function F 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 D = 1 . In contrast, when we consider the higher-dimensional small world model, each dimension D = 1 , 2 , . . . , k / 2 is considered, with the dimension ˆ D finally chosen that which provides the smallest value of F .

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 φ, d , T and λ Half 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 F 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 D = 1 ) and the higher-dimensional extension (fitted dimension D = ˆ D ). The fitted network dimension ˆ D is consistently three or less, and has median ˆ D = 2 . The table shows that allowing dimension D to exceed one allows network properties to be fit with lower rewiring rate p.

In Fig. 6 we compare fitting error (to topological properties) F and prediction error (of dynamical properties) P 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 D as well as rewiring rate p reduces fitting error F in most cases, and rarely increases error. There is a single counter-intuitive case in which fitting error F for D = ˆ D exceeds that for D = 1 . This corresponds to a situation in which statistical fluctuations cause dimension D = ˆ D > 1 to minimize F using the mean over 100 trials calculated during the fitting procedure, but this fitted dimension did not minimize F 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 P in the prediction of dynamical properties. The networks for which the original small-dimensional model provides (possibly equal) lowest F and P tend to coincide: for twenty-two networks D = 1 provides both lowest fitting error F and lowest prediction error P . 

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 R 0 .

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 personto-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 λ 1 ; and S contacting I neighbours and changing into E state with rate λ 2 . 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 γ .

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 μ = 1 / 7 , and a fourteen day convalescence period [43, 79] provides a recovery rate γ = 1 / 14 . For other parameters we consider a range of reasonable values.

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

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 is

where S, E, I and R denote the fraction of individuals in the S, E, I and R state respectively, so that S + E + I + R = 1 .

We will use the next generation matrix (NGM) [80] method to deduce the basic reproduction number. Let X and Y 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 functions

where F represents the rate of transmissions (from non-possessive to possessive compartments) and V describes transitions (between possessive or from possessive to non-possessive), such that ˙ X = F − V . 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 ( S(t ) , E(t ) , I(t ) , R (t) ) = S 0 , 0 , 0 , R 0 . 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 S 0 → 1 , R 0 → 0 at the DFE. It follows that the Jacobian matrices of F and V at the DFE are

The basic reproduction number of model (3) can be expressed in terms of spectral radius as

where ρ(M) denotes the spectral radius of the square matrix M.

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 [ X] be the number of sites that are of type X, let [ XY ] be the number of ordered X-Y pairs and, similarly, let [ XY Z] 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) are

In order to close Eqs. (4) and (5) , we must express the third-order quantities, like [ X Y Z] , [ X Y X ] , etc., in terms of the secondorder state variables [ X Y ] , [ X X ] . 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 R 0 can be estimated by introducing correlations C XY between individuals of distinct type X and Y :

This definition allows system (4) to be rewritten

By applying the NGM method as in Section 3.2.1 , we can estimate the basic reproduction number as

Assuming that exposed and infected individuals are chosen randomly at the initial time, C SE (0) ≈ 1 ≈ C SI (0) . Therefore, under PA we estimate the same basic reproduction number R 0 as under the MF method.

In this section we present the ICMC [54] 

represent the probability at time t for the node i to be S, E, I and R respectively. These probabilities satisfy the normalisation

The probabilities q E i (t) and q I i (t) for node i in state S not to be affected by any E and I neighbours are respectively

The ICMC system can therefore be expressed

We use the NGM method to deduce the basic reproduction number R 0 . Analogously to Section 3.2.1 , E i and I i , i = 1 , 2 , . . . , N are the possessive compartments, and the compartmental model can be written

, where 0 ∈ R N and i = 1 , . . . , N. Analogously to Section 3.2.1 , the disease free equilibrium (DFE) solution of ICMC method is

The Jacobian matrices of F and V at DFE point are then

where D is the identity matrix. The basic reproduction number of system (5) is

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

The fourth and final continuous state approach which we consider is the PDMC [55] method. Let X i (t) denote the probability that node i is in state X at time t, and observe that S i (t) + E i (t) + I i (t) + R i (t) = 1 . The PDMC update rules for node i

where q E i (t) , q I i (t) are given by Eq. (7) .

To identify R 0 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 i to be in states E and I are E i = E i 1 , I i = I i 1 , and S S i ≈ 1 . By Eq. (7) ,

Inserting this into system (9) , we obtain

By eliminating I i we reach

where 

In this section, we present numerical simulations illustrating the SEIR epidemic spreading model in D -dimension small world networks. Unless stated otherwise, we consider networks of size about N = 3 , 0 0 0 and mean degree k = 10 , and use epidemic parameters λ 1 = 0 . 005 , λ 2 = 0 . 01 , δ = 0 , ξ = 0 and γ = 1 / 14 , μ = 1 / 7 . Initial conditions are chosen as E(0) = 50 /N,

. 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:

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 N = 3 , 025 , dimension D = 2 , mean degree k = 10 and random rewiring rate p = 0 . 01 . Figure 8 (a,b) show the change with time of the infected and exposed node fraction I(t ) and E(t ) . These fractions are lower under MC than under other methods, and PA has the least extreme deviations from MC. Figure 8 (c,d) show f SE = [ SE] / ( k N ) and f SI = [ SI] / ( k N ) ; f XY 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 C SE , C SI for each of the five methods. The fractions of pairs of SE and SI, i.e., f [ SE] and f [ SI] , and the SE and SI correlations C [ SE] and C [ SI] 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.

We use Fig. 9 to show the final number of people affected by the epidemic, which is represented by the final affected ratio R (∞ ) = lim t→∞ R (t) , depends on infection rates λ 1 and λ 2 for a network of size N = 3 , 025 , mean degree k = 10 and dimension D = 2 , with random rewiring rate p = 0 . 01 . Under any method, the final affected ratio R (∞ ) increases as either infection rate increases. Comparing Fig. 9 (a)-(e) shows that, across the spectrum of infection rates λ 1 , λ 2 considered, the MF, ICMC and PDMC methods lead to similar final affected ratios R (∞ ) , 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 R (∞ ) varies with infection rate λ 2 for each of the five methods considered, for fixed infection rate λ 1 = 0 . 005 . 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 R 0 = 1 . 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 R (∞ ) .

In Fig. 10 we show how the final affected ratio R (∞ ) depends upon epidemic intervention parameters. Figure 10 (a) shows how R (∞ ) 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 R (∞ ) 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.

The four continuous state methods predict far stronger dependencies on epidemic parameters than are observed under MC. This is especially evident for infection rates λ 1 and λ 2 , 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 R (∞ ) , although the PA method does so least.

Now we examine the dependence of the final affected ratio R (∞ ) on parameters of the small world network model. Figure 11 presents the variation of final affected ratio R (∞ ) with mean degree k for five dimensions from D = 1 to 5. The figure reveals that mean degree k strongly affects final outcomes. Specifically, Fig. 11 (a)-(e) reveal that, in any dimension, the final affected ratio grows with mean degree k . 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 R (∞ ) compared with MC, but provides results closer to MC than are provided by the other three continuous state methods. In particular, in dimension D = 4 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 R (∞ ) grows with increasing dimension D . Thus, analyses based on MC simulations with the original, one-dimensional, small world model have the potential to underestimate the severity of an epidemic.

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 k . Instead, the discrepancy initially increases with mean degree k , reaches a peak around k = 10 , 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 D = 2 , 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 k = 10 and networks of size about N = 3 , 0 0 0 . For all five methods in dimension D = 2 , Fig. 12 (a) and (b) illustrate the variation of the final affected ratio R (∞ ) with random and distance dependent rewiring rate p respectively. Under MC we observe that the final affected ratio R (∞ ) grows with increasing rewiring rate p but, when using the continuous state methods MF, PA, ICMC and PDMC, outcomes vary little with rewiring rate p. For large random or distance rewiring rate p, the final affected ratio R (∞ ) under MC approaches that under PA but, for the other three continuous state methods, a difference persists. Figure 12 (c,d) shows how the final affected ratio R (∞ ) changes with random and distance dependent rewiring rate p under MC and in dimensions D = 1 to 5 (under the four continuous state methods, R (∞ ) varied little with rewiring rate p in any dimension). For fixed rewiring rate p < 0 . 5 , the final affected ratio R (∞ ) grows as dimension D increases from D = 1 , the dimension of the original small world model. For lower dimensions D , the final affected ratio R (∞ ) grows as rewiring rate p increases. When the random rewiring rate p is high, the final number of affected nodes varies little with dimension D , 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 ( Fig. 13 (b,c) ) or, under some conditions, slightly underestimates R (∞ ) ( Fig. 13 (a) ). In particular, the figure suggests that the discrepancies between methods are not due to finite size effects.

In Section 2 we showed that important network properties vary substantially with rewiring rate p and dimension D , 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 R (∞ ) , but the PA method does so least.

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 R (∞ ) can be characterized using network properties. To do so, we consider randomly chosen small world networks of fixed mean degree k , size between about ˜ N min = 500 and ˜ N max = 10 , 0 0 0 and rewiring rate bounded below by p min = 10 −4 ( Section 2.2 ), and compare the final affected ratio R (∞ ) X for each continuous state method X with the mean value R (∞ ) MC of final affected ratio over five MC trials. We consider 100 independently generated networks per mean degree k , and use the Kendall τ correlation coefficient to assess the statistical significance of variations of the magnitude of the deviation with changing network size N, efficiency η and capacity χ ( Section 2.1 ). In Fig. 14 we demonstrate that, for fixed mean degree k , the capacity χ governs the deviation from MC of the most accurate method, PA, more consistently than either the network size N or efficiency η. As Fig. 14 (a) shows, for fixed mean degree k = 10 , the decrease of error with capacity χ is significant at the 10 −6 level. Fig. 14 (b) does reveal that, as mean degree k 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 N and efficiency η, but only capacity χ consistently leads to correlation which is significant at the 5% level ( Fig. 14 (c) ). We have just noted that, for fixed mean degree k , the deviation of PA from MC decreases with capacity χ. However, Fig. 14 (d) shows that, conversely, when different mean degrees k are grouped, the deviation instead increases with capacity χ with Kendall τ correlation significant at the 10 −3 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. Figure 15 pertains to differences in the final affected ratios R (∞ ) 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 k = 10 , both the efficiency η and capacity χ govern deviation from MC better than the network size N for fixed mean degree k . We show in Fig. 15 (a) that, for fixed mean degree k = 10 , the decrease of error with efficiency η is significant at the 10 −13 level. As shown in Fig. 15 (a) , for mean degree varying from k = 2 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 N. However, the Kendall τ correlation is significant at the 5% level for the efficiency η for precisely the same set of mean degrees k 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 

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 di-mension 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) higherdimensional 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 R 0 . The dependence of R 0 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 R (∞ ) varies little with network size, either for MC or the four continuous state methods. Fourth, the final affected ratio R (∞ ) 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 N, but a new measure -the network capacity χ = Nη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 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 .

Dimension of spatially embedded networks

Role of fractal dimension in random walks on scale-free networks

High-dimensional Apollonian networks

Models of the small world

How to calculate the fractal dimension of a complex network: the box covering algorithm

Correlation dimension of complex networks

The correlation fractal dimension of complex networks

Random walks in varying dimensions

Random walk on random walks: higher dimensions

Synchronization of coupled oscillators on Newman-Watts small-world networks

Cluster explosive synchronization in complex networks

Small-world topology can significantly improve the performance of noisy consensus in a complex network

Diffusion and consensus in a weakly coupled network of networks

Heat flux distribution and rectification of complex networks

Influence of the degree of a complex network on heat conduction

Knowledge transmission model with consideration of self-learning mechanism in complex networks

Knowledge transmission model with differing initial transmission and retransmission process

Double transition of information spreading in a two-layered network

Review mechanism promotes knowledge transmission in complex networks

Effects of two channels on explosive information spreading

The distinct roles of initial transmission and retransmission in the persistence of knowledge in complex networks

Prisoners dilemma game model for e-commerce

Evolutionary investor sharing game on networks

Evolution of cooperation on independent networks: the influence of asymmetric information sharing updating mechanism

Synergistic effects of self-optimization and imitation rules on the evolution of cooperation in the investor sharing game

Influence of network structure on rumor propagation

A general Markov chain approach for disease and rumour spreading in complex networks

Prediction of competitive diffusion on complex networks

Steady states of epidemic spreading in small-world networks

On the stability of epidemic spreading in small-world networks: how prompt the recovery should be?

Different epidemic models on complex networks

Spreading of epidemics on scale-free networks with nonlinear infectivity

Preferential imitation can invalidate targeted subsidy policies on seasonal-influenza diseases

Precisely identifying the epidemic thresholds in real networks via asynchronous updating

Equivalence and its invalidation between non-Markovian and Markovian spreading dynamics on complex networks

Community lockdowns in social networks hardly mitigate epidemic spreading

Socio-demographic and health factors drive the epidemic progression and should guide vaccination strategies for best COVID-19 containment

Collective dynamics of 'small-world' networks

Scaling and percolation in the small-world network model

Percolation and epidemics in a two-dimensional small world

Two-dimensional small-world networks: navigation with local information

Small world and scale free model of transmission of SARS

Modelling strong control measures for epidemic propagation with networks-a COVID-19 case study

Navigation in a small world

The small-world phenomenon: an algorithmic perspective

Graph diameter in long-range percolation, Random Struct

A spatial stochastic model for rumor transmission

Dynamical systems on networks

Controllability of complex networks

Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions

Epidemic processes in complex networks

The effects of local spatial structure on epidemiological invasions

Virus spread in networks

Discrete-time Markov chain approach to contact-based disease spreading in complex networks

Epidemic dynamics in finite size scale-free networks

Accuracy of mean-field theory for dynamics on real-world networks

Non-mean-field behavior of the contact process on scale-free networks

Effect of risk perception on epidemic spreading in temporal networks

Epidemic spreading in real networks: an eigenvalue viewpoint

The second largest eigenvalue of a graph (a survey)

Spectra of Graphs

Bipartite subgraphs and the smallest eigenvalue

Max cut and the smallest eigenvalue

A characterization of the smallest eigenvalue of a graph

The largest eigenvalue of a graph: a survey

Networks: An Introduction

Chains of affection: the structure of adolescent romantic and sexual networks

Quarantine, isolation and lockdown: in context of COVID-19

Bidirectional contact tracing could dramatically improve COVID-19 control

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV-2)

Contact network models matching the dynamics of the COVID-19 spreading

A familial cluster of infection associated with the 2019 novel coronavirus indicating possible person-to-person transmission during the incubation period

Feasibility of controlling COVID-19 outbreaks by isolation of cases and contacts

Modeling the transmission dynamics of COVID-19 epidemic: a systematic review

Optimal inference of the start of COVID-19

Early dynamics of transmission and control of COVID-19: a mathematical modelling study

Spread and dynamics of the COVID-19 epidemic in Italy: effects of emergency containment measures

Mathematical modeling of the spread of the coronavirus disease 2019 (COVID-19) taking into account the undetected infections. the case of china

The construction of next-generation matrices for compartmental epidemic models

The effects of local spatial structure on epidemiological invasions

The spread of infectious diseases in spatially structured populations: an invasory pair approximation

Representing spatial interactions in simple ecological models

Reproduction numbers for epidemics on networks using pair approximation

Dynamical interplay between awareness and epidemic spreading in multiplex networks

Optimal COVID-19 quarantine and testing strategies